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parameters,  based  on  the  input  output  relations  of  the 
system  under  test,  are  developed.  The  lattice  is 
redefined  in  terms  of  the  frequency  domain  representation 
of  the  input  data.  This  proves  to  be  useful  because  it 
allows  the  input  to  be  normalized  so  that  the  lattice 
yeilds  a  consistant  set  of  parameters  independent  of  the 
test  source  characteristics.  Lastly  the  lattice  is 
redefined  in  terms  of  correlations  of  the  input  signals. 
This  results  in  a  computationally  and  storage  efficient 
lattice  algorithm. 
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ABSTRACT 

The  single  channel  autoregressive  lattice  has  been 
successfully  applied  to  problems  including  speech  analysis 
and  recognition,  spectral  analysis  and  noise  cancelling. 

More  recently  the  two  channel  autoregressive  (AR)  lattice 
has  been  exploited  for  autoregressive  moving  average  (ARMA) 
analysis  of  systems  for  modeling  and  identification. 

This  dissertation  considers  the  multichannel  AR  lattice 
when  applied  to  ARMA  systems  analysis.  Constraints  on 
lattice  parameters,  based  on  the  input  output  relations  of 
the  system  under  test,  are  developed.  The  lattice  is 
redefined  in  terms  of  the  frequency  domain  representation  of 
the  input  data.  This  proves  to  be  useful  because  it  allows 
the  input  to  be  normalized  so  that  the  lattice  yields  a 
consistant  set  of  parameters  independent  of  the  test  source 
characteristics.  Lastly  the  lattice  is  redefined  in  terms 
of  correlations  of  the  input  signals.  This  results  in  a 
computationally  and  storage  efficient  lattice  algorithm., 
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I .  INTRODUCTION 


In  many  applications  the  use  of  mathematical  models 
provides  critical  insight  into  dynamic  system  behavior  and 
is  essential  to  many  control,  simulation  and  communication 
applications.  Models  have  been  applied  in  the  areas  of 
economics,  neurophysics  and  astronomy  with  varying  degrees 

of  success.  Engineers  routinely  use  models  when  they  wish 

v 

to  predict  or  control  a  systemk  performance. 

In  this  dissertation  one  form  of  model,  the  multichannel 
lattice,  is  examined  in  detail  in  an  effort  to  gain  insight 
into  its  operation  and  to  develop  more  efficient  algorithms 
for  its  implementation  and  application  to  the  modeling  of 
multichannel  linear  and  nonlinear  autoregressive  moving 
average  (ARMA)  systems,  based  upon  a  least  mean  square  error 
approach. 

A.  THE  PROBLEM  STATEMENT 

The  output  y(k)  of  the  discrete,  linear,  time  inyariant 
system  of  Figure  1.1  may  be  described  as  the  weighted 
summation  of  the  present  and  past  values  of  the  input  u(k) 
and  past  values  of  the  output  y(k) .  This  can  be  summarised 
by  the  relation 

q  p 

v(k)  =  E  a  u(k-n)  +  E  3  y(k-n)  (1-la) 

n=0  n  n=l  n 
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defining  the  transfer  function 

Y( z)  _  a( z) 
tT(z)  ~  3  ( z ) 

where 

q  -n  P  _n 

a(z)  =  E  a  z  and  3  (z)  =  1-  E  3  z  (1-lb) 

~_n  n  n 


u(k) 


Linear  Time 
Invariant  Digital 
System 


y(k) 


Figure  1.1  General  Discrete  System 

where  a^,  0_<n_5i  and  bn,  l.fnjfP  are  "the  parameters  of  the 
system  and  k  denotes  the  sample  index. 

An  estimate  of  the  system  parameters  can  be  made  if  the 
output  signal,  y(k),  and  possibly  the  input  signal,  u(k),  of 
a  system  can  be  observed.  The  result  is  a  model  of  the 
system.  There  are  numerous  reasons  why  this  might  be  done, 
among  which  are: 

1)  The  parameters  of  the  system  are  not  known  and  some 
insight  into  the  system  function  is  desired. 
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2)  The  nominal  values  of  the  system  parameters  are 


known  and  we  wish  to  detect  a  change  in  the  system 
performance . 

Equation  (1-1)  suggests  that  a  system  output  is  predictable 
from  a  linear  combination  of  past  outputs  and  inputs. 

Because  of  this,  procedures  to  find  these  parameters  are 
often  referred  to  as  linear  prediction. 

One  way  to  construct  a  model  is  to  hypothesize  a  system 
of  the  same  form,  namely 

A  /N  P  /V  /N 

y(k)  =  I  a  u(k-n)  +  I  b  y(k-n)  (1-2) 

n=0  n  n=l  n 

A  A 

and  adjust  the  model  parameters  a  and  b  so  that  the 

n  n 

difference  between  y(k)  and  y(k)  is  minimized  according  to  some 
criteria.  Figure  1.2  shows  this  method.  The  symbol  is 

used  to  indicate  an  estimated  value.  From  Figure  (1.2)  and 
(1-1)  and  (1-2) 

e(k)  =  y(k)  -  y(k) 

q  P  k 

-  Z  (a  -a  )u(k-n)  +  Z  3  y(k-n)  -Zb  y(k-n) 
n-G  n=l  n=l 

A 

Since  y(k)  is  itself  dependent  on  model  parameters,  and 
a  minimum  mean  square  output  error  criterion  results  in  a 
set  of  nonlinear  equations,  the  solution  of  this  equation 

A 

for  model  parameters  is  generally  intractable  when  b  , 

n 

l£n<p  is  not  zero. 
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Figure  1.2  Direct  System  Modeling 

Because  of  this,  other  analysis  procedures  have  been 

developed.  A  very  powerful  model  is  the  equation  error 

formulation,  which  entails  two  models;  the  analysis  model, 

with  parameters  a  and  b  which  are  estimates  of  the  svstem 

n  n 

parameters,  and  a  synthesis  model  using  a^  and  bn  to  emulate 
the  system. 

The  output  of  the  equation  error  analysis  model  is 
written  as 

q  p 

y(k)  =  E  a  u(k-n)  +  E  b  y(k-n)  (1-3) 

_ n  n  .  n 

n=0  n=l 
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where  the  estimated  output  is  now  a  function  of  the  system 
past  and  present  input  and  system  past  outputs  (instead  of 
the  model  outputs).  This  form  is  shown  as  Figure  1.3.  With 
the  equation  error  approach  a  MMSE  solution  for  the  model 
parameters  is  a  set  of  simultaneous  linear  equations  with  a 
unique  minimum.  The  synthesis  model  for  the  equation  error 
approach  is  of  the  form  of  (1-2),  using  the  estimated  system 
parameters  found  from  the  analysis  model.  The  equation 
error  of  Figure  1.3  is  given  by 

-  q  P 

e(k)  =  y(k)  -  y(k)  =  E  ( a  -a  )u(k-n)  +  E  (3  -b  )y(k-n) 

n  =  0  n  n  n=l  n  n 

(1-4) 

which  is  linear  in  the  model  parameters.  When  a  =a  and 

n  n 

6n=bn,  e(k)=0.  In  terms  of  the  z-transform, ( 1-4 )  can  be 
written  as 

E(z)  = 


where 

A(z)  = 


-A(z)U(z)  +  3(z)Y(z) 


C-A(z)  +  |[|y  cx(z>]  U(z) 


(1-5) 


q  p 

E  a  z_n  and  3(z)  =1-  E  b  z  n. 
«  n  _  t  n 

n=0  n=l 
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u(k) 


System  To  Be 
Modeled 


y(k) 


Figure  1.3  Equation  Error  Analysis  Model 

Three  cases  are  examined  in  this  dissertation. 

1)  The  all  pole,  or  autoregressive  (AR) ,  model 
an  =  °>  i_5*iq 

2)  The  all  zero,  or  moving  average  (MA) ,  model 
bn  =  0 ,  l_<n_<p 

3)  The  pole-zero,  or  autoregressive  moving  average 
(ARMA),  model 
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an  t  0,  for  some  n,  l£n_<q 
i  0 ,  for  some  n,  l<n<p 

The  objective  will  be  to  develop  procedures  and  algo¬ 
rithms  that  allow  more  efficient  calculation  of  these  models 
for  application  to  system  performance  evaluation  and 
identification. 

B.  OVERVIEW  OF  THE  DISSERTATION 

The  following  chapters  are  first  a  review  of  some 
contemporary  thoughts  on  linear  systems  modeling  using 
autoregressive,  moving  average  and  autoregressive  moving 
average  modeling,  with  their  lattice  formulation,  followed 
by  some  new  lattice  configurations  centered  upon  the  ARMA 
lattice  configurations. 

Yule  [19]  first  developed  the  autoregressive  model  in  a 
paper  on  sunspot  analysis.  Because  of  its  wide  ranging 
applications,  autoregressive  modeling  has  received  extensive 
examination  since.  Chapter  II. A  begins  by  reviewing  the 
relevant  theory  of  the  least  mean  square  solution  of  the  AR 
equations.  This  solution  results  in  a  set  of  simultaneous 
equations  that  are  the  equivalent  of  the  Yule-Walker  [9] 
equations.  Their  resolution  requires  the  solution  of  the 
matrix  equation  Rh  =  r,  where  R  is  a  symmetric  and  Toeplitz 
correlation  matrix.  A  number  of  methods  have  been  proposed 
for  the  solution  of  these  equations.  Matrix  inversion 
methods  can  be  used,  but  their  computational  complexity  is  of 
concern  for  higher  order  problems.  Adaptive  schemes  [18] 


u. ; 
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have  been  developed  that  update  an  estimate  of  the  solution 
as  each  data  point  is  taken  into  account.  Levinson  [6] 
presented  an  elegant  and  efficient  recursive  method  of 
performing  this  inversion  by  finding  an  (n+l)-th  order 
solution  as  an  adjustment  of  the  n-th  order  solution.  Durbin 
[2]  later  recognized  that  the  elements  of  the  vector  r  were 
also  elements  of  the  matrix  R.  This  was  used  to  improve 
the  computational  efficiency  of  Levinson's  algorithm.  This 
recursive  solution  can  be  implemented  as  a  lattice  structure 
that  has  been  found  to  offer  an  accurate  and  efficient 
solution  to  the  AR  modeling  problem.  More  recently  Burg 
Cl]  has  shown  that  this  lattice  is  a  maximum  entropy 
solution,  proving  that  the  lattice  solution  is  optimum. 

In  Section  II. B  the  moving  average  model  is  similarly 
discussed.  The  result  of  a  least  mean  square  solution  is 
the  digital  form  of  the  Weiner-Hopf  [17]  equations  and  also 
requires  the  inversion  of  a  symmetric  and  Toeplitz  matrix. 
Perry  [16]  has  shown  that  this  can  be  efficiently  solved 
using  Levinson's  recursion  and  can  also  be  restructured  to 
form  a  lattice  implementation,  part  of  which  is  the  AR  model. 

Sections  II. A  and  II. B  deal  with  single  channel  (single 
input,  single  output)  models.  Section  II. C  discusses 
multiple  input,  multiple  output  generalizations  of  the  AR 
and  MA  equations  and  structures  that  Robinson  [16]  has 
shown  can  be  solved  as  a  multichannel  generalization  of  the 
single  channel  structures. 
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Chapter  III  is  a  discussion  of  ARMA  modeling.  Unlike 
the  AR  and  MA  models,  the  ARMA  model  represents  both  poles 
and  zeros  and  so  will  offer  a  lower  order  solution  than 
either  the  AR  or  MA  models  for  zero-pole  systems.  In  III. A 
the  least  mean  square  solution  of  the  ARMA  problem  is 
developed.  A  direct  application  of  the  ARMA  equations 
results  in  a  set  of  nonlinear,  multimodal  equations.  This 
problem  was  greatly  simplified  when  the  equation  error 
formulation  [5,  12]  was  found  to  result  in  linear  equations. 
Then  in  Section  III.B  it  is  shown  how  Perry  [16]  recon¬ 
figured,  with  appropriate  assumptions,  the  ARMA  problem 
so  that  it  could  be  solved  using  a  structure  identical  to 
the  two  channel  AR  lattice.  This  structure  offers  the 
advantages  of  a  recursive  in  order  solution  and  results 
in  more  accurate  solutions  than  batch  methods  that  use 
more  direct  matrix  inversions . 

Lattices  have,  in  the  past,  used  data  in  the  form  of  a 
discrete  time  series  as  a  block  of  data,  as  discussed  in 
Chapter  II  and  III. A  and  III.B.  Griffiths  [3]  and  later 
Morf  [13]  developed  adaptive  algorithms  that  still  use 
discrete  time  series  data  but  update  an  estimate  of  the 
lattice  parameters  as  each  data  point  enters  the  lattice 
structure.  In  Section  III.C  a  new  formulation  for  deter¬ 
mining  the  lattice  parameters  is  presented.  This  approach, 
although  fundamental,  has  generally  not  been  exploited, 
particularly  for  the  multichannel  case.  In  this  structure 
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the  data  presented  to  the  lattice  is  the  frequency  domain 
spectrum  of  the  signal  to  be  analysed.  The  power  spectral 
desnity  of  the  error  spectrum  is  minimized  and  the  resultant 
lattice  is  equivalent  to  that  used  with  time  domain  data. 

This  new  lattice  structure  is  useful  because  it  allows  the 
spectram  content  of  the  input  signals  to  be  observed  and 
presents  the  possibility  of  adjusting  the  input  spectra 
to  compensate  for  a  nonwhite  driving  source.  This  configuration 
also  permits  the  lattice  to  be  used  for  filter  synthesis. 

Section  III.D  discusses  the  lattice  ARMA  model  when  applied 
to  system  identification.  The  ARMA  lattice  identifier  can  be 
used  for  system  performance  evaluation.  In  this  application 
the  model  parameters  obtained  as  a  result  of  a  test  are  com¬ 
pared  with  model  parameters  obtained  from  a  system  functioning 
normally  and  functioning  with  known  faults.  This  comparison 
can  reveal  if  the  system  under  test  is  operating  satisfactorily 
and  possible  identify  a  class  of  faults  within  the  system. 

The  lattice  parameters  fully  identify,  in  a  least  mean  square 
sense,  both  the  input  and  output  of  the  system  to  be  modeled. 

For  modeling  it  is  generally  desirable  that  the  input  signal 
be  white  (have  a  flat  spectrum)  so  that  the  signal  energy  is 
equally  distributed  across  the  spectrum  of  interest.  Some  new 
relations  among  the  lattice  parameters  are  developed,  which, 
with  some  previously  known  relationships,  reduce  the  number 
of  parameters  needed  by  one  half  when  the  input  is  properly 
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conditioned.  Some  new  methods  are  then  derived  that  allow 
preemphasis  or  normalization  of  the  system  input  to  correct 
for  a  nonwhite  input  signal.  This  is  important  because 
lattice  parameters,  while  modeling  the  system,  are  still  a 
function  of  the  input  signals  second  order  statistics.  This 
is  useful  even  when  the  test  source  is  purportedly  white, 
since  the  lattice  must  function  with  estimates  of  signal 
parameters  obtained  from  finite  length  data  signals. 

Section  III.F  introduces  another  application  of  the  new 
frequency  domain  lattice.  The  lattice  can  be  used  to  experi¬ 
mentally  synthesize  a  filter  to  meet  a  desired  specification 
by  providing  the  desired  response  as  the  input  of  the 
lattice.  The  model  generated  is  a  minimum  mean  square  error 
approximation  of  the  desired  specification. 

In  Chapter  IV  and  Appendix  A  a  new  approach  to  both  the 
single  channel  and  two  channel  lattice  is  developed.  This 
method  recognizes  that  the  error  signals  of  the  lattice  are 
only  linear  combinations  of  the  input  signals  and  delayed 
versions  of  the  input  signals.  From  this  observation  the 
equations  for  the  lattice  parameters  are  restructured  so  that 
only  linear  combinations  of  the  input  correlation  functions 
are  needed  to  solve  the  lattice.  This  is  more  computational 
and  storage  efficient  than  other  lattice  methods  since  it  is 
not  now  necessary  to  generate  updated  error  signals  at  each 
stage.  It  retains  the  advantages  of  the  lattice  of  orthogonality 
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and  maximum  entropy.  It  offers  the  new  advantage  that  data 
need  not  be  artifically  windowed  to  force  the  correlation 
matrix  to  be  Toeplitz,  as  is  required  by  Levinson's 
algorithm.  This  may  lead,  with  some  types  of  data,  to  a 
more  accurate  solution.  This  method  is  extendable  to  higher 
dimensioned  lattices  as  are  used  for  nonlinear  and  multiple 
input  multiple  output  system  modeling. 

Chapter  V  presents  a  summary  of  the  new  results  in  this 
dissertation.  Included  is  a  discussion  of  the  application 
and  implementation  of  these  new  concepts,  along  with  questions 
for  future  research. 
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II.  LINEAR  SYSTEM  MODELING  AND  IDENTIFICATION 

The  autoregressive,  moving  average,  and  autoregressive 
moving  average  models  have  all  been  successfully  used  for 
system  modeling  and  analysis.  The  AR  model  is  the  model  of 
choice  when  only  system  output  is  available.  The  MA  model 
can  give  lower  MSE  if  an  input  signal  is  available,  but 
cannot  generate  an  infinite  impulse  response.  These  models 
will  be  examined  in  greater  detail  and  computationally 
efficient  algorithms  for  their  solution  will  be  developed. 
Both  single  input  and  single  output  and  multiple  channel 
implementations  are  discussed. 

A.  THE  AUTOREGRESSIVE  MODEL 

For  the  autoregressive  model  the  assumption  is  made  that 
the  present  output  of  a  system  can  be  estimated  from  a 
weighted  summation  of  the  past  values  of  the  output.  This  can 
be  expressed  as 
P 

y(k)  =  I  b  y(k-n)  (2-1) 

n=l 

or  in  matrix  form 

y (k)  =  bT£(k)  (2-2) 
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The  mean  squared  prediction  error  as  a  function  of  the 
weights  bn  is  found  to  be 

E  ■  feTSyy&  -  2  -^-yy  +  V0’  <2-4> 

where  in  general  R  (n)  =  t  {x(k)w(k+n)} ,  r  =  e {x(k)w(k+n)>  , 

xw  — xw  — 

T 

R^w  =  e{x(k)w  (k)}  and  e{*}  designates  expectation. 

In  this  case, 

V<0)  •••  V('p+1) 

R"  (p-1)  ...'R  (0) 

_yy  yy  J 

Eyy  =  [V(1)  •••  Ryy(p)] 


The  solution  of  (2-4)  for  the  optimum  weights  is  given  by 


b  . 
-opt 


(2-5) 


which  are  commonly  called  the  normal  equations.  If  this  is 
substituted  into  (2-4)  then  the  minimum  error  is  found  to  be 


^min  ”  ^yy^^  “  ^opt  — yy 


(2-6) 


(2-3)  may  be  written  in  the  z  domain  as 


E(z) 

TTZ)  =  1 


P 

E 

n=l 


bn2 


-n 


B(z) 


(2-7) 


which  is  the  relationship  of  error  signal  output  to  the 
input  signal. 


22 


Consider  the  all  pole  filter  of  Figure  2.2  with  the 
input  signal  u(k)  and  output  y(k).  The  relationship 
between  input  and  output  can  be  written  in  the  z  domain 


Y(z) 

OTzT 


H(z)  = 


1-  I  b_z 


n=l  “ 

or  in  the  time  domain 

P 

y(k)  =  anu(k)  +  I  b  y(k-n) 
0  n=l  n 


( 2-8a) 


( 2  —  8b ) 


u(k)  a0/^\ 


y(k) 

1 - »z_1— , 

Figure  2.2  All  Pole  Filter 

If  the  output  y(k)  of  the  all  pole  filter  is  applied  to  the 
input  of  the  prediction  error  model  of  Figure  2.1  the  output 
error  will  be  e(k)  =  aQu(k).  The  transfer  function  of  this 
cascade  arrangement  from  the  input  of  the  all  pole  filter  to 
the  error  output  of  the  predictor  will  be  aQ.  This  shows  that 
if  u(k)  is  white  noise  the  effect  of  the  prediction  error 
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model  is,  except  for  a  gain  term  a„  ,  to  reverse  the  action 
of  the  all  pole  filter  or  to  whiten  it  *  s  output.  For  this 
reason  the  prediction  error  model  is  often  referred  to  as  a 
whitening  filter  or  an  inverse  filter,  and  makes  it  useful 
as  a  spectrum  analyzer.  It  should  be  noted  that  a^  establishes 
a  lower  bound  on  the  MSE  of  the  AR  predictor.  If  u(k)  is  a 
unit  impulse,  u(k)  =  5  (k) ,  the  model  may  replicate  the 
response  y(k)  perfectly,  but  can  never  predict  the  arrival 
of  the  input  a^6 (k)  since  the  predictor  is  causal,  and  there 
is  always  this  minimal  error. 

The  relationship  between  an  all  pole  filter  of  Figure 
2.2  and  the  analysis  model  of  (2-1)  suggests  that  the  all 
pole  filter  can  be  used  as  a  synthesis  model  for  generating 
the  signal  y(k)  using  aQu(k)  as  the  input  signal  with  the 
coefficients  obtained  from  the  prediction  error  analysis 
model.  It  is  only  necessary  to  find  a  suitable  gain  term 
aQ.  Equation  (2-3)  can  be  rewritten 


y(k)  =  e(k)  +  b"^(k)  (2-9) 

Comparing  this  with  ( 2  —  8 b )  it  is  seen  that  aQu(k)  =  e(k). 
3y  requiring  that  the  total  energy  oi  the  output  from  the 
linear  predictor  be  equal  to  the  total  energy  of  the 
input  to  the  all  pole  filter,  a^  is  specified  by 


2  _  z  {e2  (k)} 
0  '  Ruu(0) 


(2-10) 
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This  leads  to  a  2  domain  transfer  function  representation 
for  the  synthesis  model  of 


H(z) 


ao 

3Tz7 


(2-11) 


The  impulse  response  of  (2-11)  can  be  of  infinite  dura¬ 
tion.  This  means  that  a  relatively  low  order  AR  model  may 
emulate  an  infinite  impulse  response  system  to  an 
acceptable  accuracy. 

In  the  problem  as  stated  so  far  it  has  been  tacitly 

assumed  that  the  desired  order  for  the  model  was  somehow 

known.  In  practice  this  is  often  not  the  case.  Instead, 

(2-5)  is  solved  for  some  order  and  the  MSE  determined.  If 

this  MSE  is  unacceptable,  then  (2-5)  is  again  solved  for  a 

higher  order.  Solution  of  the  simultaneous  equations  of 

3  2 

(2-5)  requires,  by  Gaussian  elimination,  p  /3  +  0(p“) 
operations,  so  that  a  repetitive  solution  of  the  normal 
equations  can  become  computationally  ponderous.  There  is, 
however,  a  high  degree  a  symmetry  in  these  equations  and 
the  exploitation  of  this  symmetry  can  give  a  more  efficient 
solution  of  (2-5). 

1 .  A  Recursive  In  Order  Solution 

The  normal  equations  (2-5)  for  the  optimum  weights 
can  be  written  in  matrix  form 
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R  (0)  R  (-1)  . . .  R,  (-p+l> 

yy  yy  yy 

bi~ 

R  (1) 
yy 

R  (1)  R  (0) 
yy  yy 

•  •  • 

•  •  • 

b2 

V2) 

«  •  • 

R  (p-1)  ...  R  (0) 

Lyy  yy 

_bp  _ 

_Vp). 

The  autocorrelation  matrix  is  Toeplitz  (the  elements  along 
any  diagonal  are  equal)  if  the  signals  are  stationary  which 
is  a  basic  assumption  of  this  work. 

Calculation  of  the  autocorrelation  functions  now 
deserves  some  comment.  Computation  of  the  exact  discrete 
correlation  requires  a  summation  of  infinite  duration 
signals . 

CO 

R  (N)  =  l  y(k)y(k+N) 

yy  k=-°° 


The  result  is  an  even  function,  or 


R  (N) 

yy 


R  ( -N) 

yy 


so  that  the  autocorrelation  matrix  is  symmetric  and 
Toeplitz.  In  the  absence  of  an  infinite  duration  signal, 
the  more  common  case,  an  estimate  of  the  correlation  func¬ 
tion  must  be  made.  If  a  window  function  is  used  and  the 
data  is  taken  to  be  zero  outside  the  window  then  the 
estimated  correlation  matrix  will  be  symmetric  when  the  data 
is  averaged  over  the  width  of  the  window.  If  the  data  is 
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not  taken  to  be  zero  outside  the  window  and  averaging  is 
done  over  a  finite  number  of  points,  the  estimated 
correlation  matrix  is  in  general  not  symmetric  unless  a 
large  number  of  data  points  are  taken  because  of  end 
effects.  In  the  derivations  which  follow  symmetric  Toeplitz 
matrices  are  assumed.  It  is  noted  that  the  final  algorithms 
can  be  used  without  assuming  symmetric  Toeplitz  correlation 
matrices,  and  apparently  work  just  as  well.  In  fact  it  is 
noted  that  Griffith’s  time  adaptive  technique  [3]  is 
nonstationary  and  yet  appears  to  work  very  well.  This  has 
not  been  explained  in  the  literature.  If  the  observed  signal 
is  ergodic  an  estimation  of  the  correlation  function  can  be 
made  using  a  finite  interval  of  data. 

When  the  data  is  windowed  so  that  the  autocorrelation 
matrix  is  Toeplitz  and  symmetric,  the  method  used  to  solve  the 
resultant  normal  equations  is  often  referred  to  in  the 
literature  [9]  as  an  autocorrelation  method  of  solution. 

When  the  data  is  not  so  windowed  and  the  matrix  is  not 
Toeplitz,  the  problem  solution  is  said  to  be  a  covariance 
method.  Either  method  results  in  an  estimate  of  the  true 
values  of  the  correlation  matrix,  and  it  is  generally 
unclear  that  one  method  will  be  more  accurate  than  the  other. 

The  Levinson  algorithm  relates  the  (n+l)-th  order 
solution  of  a  set  of  equations  of  the  form  of  (2-12)  to  the 
n-th  order  solution.  The  assumption  is  made  that  b^n+'*'\ 


27 


where  the  superscript  denotes  the  order  of  the  solution,  can 
be  found  from  by  adding  a  correction  to  each  term  and 

one  new  term,  or 


b(n+1)  =  -  + 


.(n+1) 


where  b 


(2-13) 


To  find  b^n+^  parameters  and  k^n+^  must  be  computed. 

Permuted  versions  of  the  vectors  b^n^  and  r^n^  may  be 

-yy 

defined  by  reversing  the  order  of  their  elements. 


„<n) 

-yy 


R  (1) 

_yy  . 


=  ;(n) 
-yy 


(2-14) 


The  symmetry  of  the  autocorrelation  matrix  allows  (2-12)  with 
p=n  namely  Ryy^b^n^  =  r^^  ,  to  be  rewritten  using  (2-14)  as 


n(n)r;(n)  _  (n) 

Kb  -  P 

-yy  -  -yy 


(2-15) 


The  matrix  R^n+^  is  recognized  as  the  matrix  R^n^  with  one 

-yy  -yy 

new  row  and  column  and  may  be  written  as 


r  r  (-nfj  r  1 

R<n)  j r  (-D  -yy  i  o 

n(n+i )  -yy  ,  ^  . — 

-yy  '  .  >  "  (  xT1 

R  (n) . . . R  ( 1 ) i R  (0)  pkn;  I  R  (0) 

L  yy  yy  yy  _  -yy  1  yy 
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since  the  data  has  been  properly  windowed  to  ensure  that  the 
autocorrelation  function  estimate  is  an  even  function  of 
time.  The  foregoing  and  (2-13)  may  be  substituted  into  an 
(n+l)-th  order  version  of  (2-12)  to  give 


R(n) 

-yy 

1  (n) 

j  -yy 

.(n).-(nf 
b  H 

1^  1 

<<  /-v 

^  o 

_ 1 

p‘->T 

-yy 

|  R  (0) 

1  yy 

k(n*l) 

J 

R  (n+1) 

yy 

Multiplying  the  matrices  and  substituting  (2-15)  yields 


r(n)  _  _^(n)^(n+l) 


(2-16a) 


r^*i)  -  4;>T  ^(n> 

R( 0 )  -  b(n)T  r(n) 

yy  -  -yy 

Substituting  (2-16a)  into  (2-13)  yields 


b(n+1)  = 

b(n)" 

+  ktn+i> 

'b(n)~ 

1 

o 

i 

l _ 

where 

k(n+1)  = 


(2-16b) 


(2-17) 


To  find  b^n+1^  from  b^n^  it  is  only  necessary  to 
determine  ,  which  can  be  found  from  (2-16b)  which  uses 

R  (n+1)  and  other  terms  from  the  previous  solution.  It  is 

yy 

noted  by  comparison  with  (2-6)  that  the  denominator  of  (2-17) 
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is  the  MMSE  of  the  (n)-th  order  solution  for  the  normal 
equations.  Equations  (2-13),  (2-14),  (2-16)  and  (2-17)  can 
be  applied  recursively,  starting  with  order  1,  until  a  model 
results  with  the  desired  accuracy. 

2 .  An  AR  Lattice  Structure 

The  vector  b^n^  can  be  written  in  the  transform 

domain 


B(n)(z) 


1  - 


n 

E 

i  =  l 


,(n)  -i 
bi  z 


=  1  - 


h(n)  -1 
bi  2 


-b(n)z'n 

n 


(2-18) 


and  the  reversed  in  order  version  b 


(n) 


can  be  written 


B(n)(z)  = 


z-nB(n)(z*1) 


(2-19) 


It  should  be  recalled  that  B(z)  is  the  transfer  function  of 
the  prediction  error  model.  The  overbar  here  is  used  to 
indicate  the  reversal  in  order,  and  will  be  used  throughout 
this  dissertation  to  signify  a  backward  signal,  or  a  function 
or  martix  associated  with  a  backward  signal. 

A  recursive  solution  for  the  transfer  function  3(z) 
can  be  found  by  substituting  (2-17)  into  an  (n+1)  order 
version  of  (2-18) 

„  ( n+1) ,  ,  _  ,  r  —1—2  -n-1-,.  (n+1) 

3  ( z)  =  1  -  Lz  z  ...z  Jb 


30 


The  forward  prediction  error  is  the  difference  in  the  signal 
and  a  predicted  signal  derived  from  a  weighted  summation  of 
the  past  n  samples  of  the  signal.  The  backward  prediction 
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error  is  the  difference  in  the  signal  at  time  (k-n)  and 
its  predicted  value  based  on  a  weighted  summation  of  the  n 
samples  to  follow.  Figure  2.3  is  an  illustration  of  this 
forward  and  backward  prediction. 


Figure  2.3  Illustration  of  Forward  and  Backward 
Prediction  Using  n  Points 


By  multiplying  (2-20)  and  (2-21)  by  Y(z),  E(z)  and 
E(z)  may  also  be  found  recursively. 

E(n+1)(z)  =  E(n)(z)  -  k(n+1)z'1  E(n)(z)  (2-23) 


ZCn^\z)  =  z’1  £(n)(z)  -  k(n+1)E(n)Cz) 


( 2-23a) 


In  the  time  domain  these  two  equations  become 


e(n+1)(k)  =  e(n)(k)  -  k(n+1)e(n) (k-1)  (2-24) 

— (n+i) (k)  _  e(n)(k-I)  -  k(n+1)e(n)(k)  (2-25) 

/v 

Since  the  zero  order  prediction  of  a  signal  is  y(k)  =  0  (no 
prediction  at  all)  the  zero  order  error  is 

e(0)(k)  =  e(0)(k)  =  y(k)  (2-26) 

and  equations  (2-24),  (2-25)  and  (2-26)  can  be  implemented  as 
a  lattice  structure  as  is  shown  in  Figure  2.4. 


Figure  2.4  Second  Order  AR  Lattice  Model 


The  lattice  parameter  k 


can  be  found  from  (2-16b) 


(as  was  done  for  the  Levinson  algorithm)  or  can  be  found  by 
minimizing  the  expected  value  of  the  squared  error  of  (2-24) 
or  (2-25).  Minimizing  square  of  the  the  expected  value  of 
(2-24)  with  respect  to  k  yields 


k 


(n+1) 

F 


e{e(n)(k)e(n)(k-l)} 
e  {e(n)(k-l)e(n)(k-D} 


(2-27) 


and  similarly  for  (2-25) 


(n+1) 

B 


e{e(n)(k)i(n)(k-l)} 

e{e(n)(k)e(n)(k)} 


(2-28) 


These  alternative  solutions  are  referred  to  as  the 
forward  and  backward  methods  and  can  be  shown  [8]  to  be 
equivalent  to  (2-16b).  From  (2-20)  and  (2-21)  it  is  seen 
that  the  forward  and  backward  prediction  transfer  functions 
are  equal  so  (2-27)  and  (2-28)  must  be  equivalent.  However, 
when  implemented,  the  expectations  will  be  estimated  from 
time  averages  and  will  generally  not  be  equal.  This  has 
led  to  two  alternative  solutions,  the  first  a  harmonic 
mean  of  (2-20)  and  (2-21)  due  to  Burg  [1] 


,(n+l) 

"BURG 


£{e(n)(k)e(n)(k-l)} 


j£{e(n)  (k)+e(n)  (k-1) } 


(2-29) 


and  another  from  Itakaura  and  Saito  [4]  which  is  the  geometric 
mean  of  (2-27)  and  (2-28) 


(n+1)  e{e(n)(k)e(n)(k-l)} 


e{e(n)  (k)}e{I(n)  (k-1)} 


(2-30) 


Stability  of  the  result  is  ensured  if  all  the  poles  of  the 
transform  domain  transfer  function  lie  within  the  unit 
circle.  Markel  and  Gray  [11}  have  shown  that  the  necessary 
and  sufficient  condition  for  this  is 


(2-31) 


Since  (2-30)  is  of  the  form  of  a  normalized  correlation  it 
is  of  necessity  always  bounded  by  unity,  and  since  the 
magnitude  of  a  geometric  mean  is  an  upper  bound  of  the 
magnitude  of  a  harmonic  mean,  (2-30)  and  (2-29)  are 
guaranteed  to  result  in  stable  solutions.  There  is  no  such 
guarantee  with  the  forward  (2-27)  and  backward  (2-28)  methods, 
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A  most  important  property  of  the  lattice  structure 
is  the  orthogonality  of  the  backward  error  signals, 
described  by 


e  {e(i)(k)i(  j)(k)> 


(2-32) 


Indeed,  the  lattice  structure  can  be  derived  [11]  starting 
with  this  property.  This  orthogonality  decouples  successive 
stages  of  the  lattice,  ensuring  that  the  addition  of  stages 
in  the  recursive  solution  of  the  filter  will  not  affect 
preceding  stages,  and  permitting  the  derivation  of  (2-27) 
and  (2-28).  Furthermore,  [8]  orthogonality  allows  for  a 
recursive  calculation  of  the  error  energy  at  successive 
stages  by 

e{e(n+1)(k)  e(n+1)(k)}  =  e{e(n+1)(k)  e(n+1)(k)}  =  p(n+1) 

?(n+l)  _  p(n)(1_k(n+l)2)  (2-33) 

Noting  that 

P<0)  =  e(e(0)(k)  e(0)(k)}  =  Ryy(0)  (2-34) 

(2-33)  and  (2-34)  may  be  used  for  the  denominator  of  either 
( 2-16b)  ,  (2-27)  or  (2-28). 
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Once  (2-17)  is  solved,  the  vector  b^n^  may  be 
implemented  as  the  all  pole  filter  of  Figure  2.2,  or  other 
equivalent  arrangements,  to  produce  a  synthesis  model. 

The  lattice  synthesis  structure  may  be  generated  by  (2-25)  and 
rearranging  (2-24). 

e(n)(k)  =  e(n+1)(k)  +  k(n+1)  £(n)(k-l)  (2-35) 

This  equation  can  be  realized  by  the  structure  of  Figure  2.5 
and  implements  the  transfer  function  1/B^n^(z). 


Figure  2.5  Second  Order  AR  Lattice  Synthesis  Structure 


B.  THE  MOVING  AVERAGE  MODEL 

The  moving  average  model  is  based  on  the  assumption  that 
the  present  value  of  a  system  output  can  be  estimated  by  a 
weighted  summation  of  the  present  and  past  q  values  of  the 
system  input.  In  equation  form  this  is 
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(2-36) 


~  H  T 

y(k)  =  Z  a(n)  u(k-n)  =  a  u(k) 
n  =  0 

where 

a?  =  [a(0)  a(l)  ...  a(q)] 

uT(k)=  Cu(k)  u(k-l)  ...  u(k-q)  ]  (2-37) 

It  should  be  noted  that  the  vectors  of  (2-37)  are  (q+1) 
elements  long. 

Defining  the  error  of  the  model  as  the  difference  in  the 
estimated  value  and  the  actual  value  of  the  system  output, 
or 


e(k)  =  y(k)  -  y(k) 


(2-38) 


the  mean  squared  error  E  is  found  to  be 

E  =  a1 R  a  -  2  aT  r  +  R  (0) 

—  -uu—  -  — uy  yy 


(2-39) 


where  R  is  (q+l)x(q+l)  and  r  is  a  (a+1)  element  vector. 
When  this  is  minimized  with  respect  to  the  coefficients  of  a 
the  result  is  the  normal  equations  of  the  moving  average 
problem 


R  a  = 
— uu— 


r 

— uy 


(2-40) 
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From  this  the  optimum  value  for  a  is 


^opt 


(2-41) 


while  the  MMSE  is  found  by  substituting  (2-41)  into  (2-39) 


E  .  =  R  (0)  -  a  .  r 

min  yy  —opt  — uy 


(2-42) 


Figure  2.6  Moving  Average  Modeling 


Figure  2.6  shows  the  block  diagram  of  the  implementation 
of  the  MA  model.  The  transfer  function  of  the  model  is,  in 
the  z  domain 

q  -i 

A(z)  =  £  a(n)z  (2-43) 

n  =  0 
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and  is  often  referred  to  as  an  all  zero  model.  As  such  its 
impulse  response  is  finite  and  in  all  cases  the  model  is 
stable . 

1 .  A  Recursive  In  Order  Solution 

Equation  (2-40)  is  similar  in  structure  to  (2-12). 

R  is  Toeplitz  as  is  R  .  of  (2-12),  but  the  elements  of  r 
— uu  — yy  — uy 

of  (2-40)  are  not  also  elements  of  the  matrix  R  .  Still 

-uu 

Levinson's  algorithm  may  be  applied. 

Equation  (2-40)  is  written  for  order  n  as 


R(n)  a(n)  =  r(n) 
— uu  —  — uy 


(2-44) 


and  the  assumption  is  made  that  a/n+1^  is  a  function  of 
and  a  correction  vector. 


(n+1) 


a(n)" 

+ 

y (n+l) 

1  O  1 
1 
1 

1  J 

g(n*l) 

(2-45) 


Define  vectors 


r<n+1)  = 

— uu 

fR  (1) 
uu 

• 

„(n  +  l ) 

’  £uu 

fR  (n+lfl 
uu 

R  (n+1) 
_uu  J 

R  (1) 

_  uu 

-  p(n+1> 

— uu 


(2-46) 
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and  expand  an  (n+1)  order  version  of  (2-44)  with  (2-45) 
substituted. 


R(n) 

— uu 

i  p<"+in 

'  AiU 

.  |  _  _ 

a(n)  +  x(n+1)‘ 

'r(n)  1 

-uy 

+  - 

r 

(n+l)T 

p 

1 

i  R  (0) 

L_s<ntl) 

R(n+1) 

Lruu 

1  uu 

L  uv  A 

When  this  is  multiplied  two  equations  result 

R(n)a(n)  +  R(n)  (n+1)  +  (n+1)  (n+1)  _  p(n) 

— uu  —  — uu  —  —  uu  °  ~  — uy 


(ntl)T  (n)  +  (n.X)T  (n+l>  +  R 
— UU  -  -uu  -  uu  6 


(n+1) 


Ruy(n+1) 


(2-48) 


The  first  and  last  terms  of  (2-47)  are  equal  so  both  may  be 
cancelled  leaving,  when  solved  for 


(n+1)  _  p(n) 

Y  -  -K 

—  — uu 


-1 


(n+1)  (n+1) 

-uu  s 


(2-49) 


Finding  R^j^ 


-1 


is  comparable  to  (2-5),  the  (n+l)-st  order 


autoregression  of  the  input  signal  u(k).  So  to  derive  the 


MA  model,  the  AR  model  of  the  input  must  first  be  solved. 
When  this  is  done  a  new  vector  f^n+^  can  be  defined  where 
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f Cn+1) 


w(n)'1  (n+1) 

R  p 

— UU  MiU 


’  ( n+ 1 )~ 
n+1 


(n+1) 


(2-50) 


is  the  result  of  the  autocorrelation  analysis.  Now  is 

defined  as 


(n+1)  ^(n+1)  (n+1) 

X  =  -  £  g 


(2-51) 


and  this  can  be  substituted  into  (2-48)  and  solved  for 


g 


(n+1) 


g<n+1) 


R 


(n+1) 
XLX _ 


R  (0) 
uu 


(n+l)T  (n) 

P  3 

— uu _ — _ 

p(n+l)T  f (n+1) 


(2-52) 


2 .  A  MA  Lattice  Structure 

The  Levinson  solution  of  (2-49)  can  be  restructured 
to  give  a  lattice  configuration.  The  transfer  function  of 
the  MA  model  (2-43)  can  be  written  recursively  in  the  z 
domain  using  (2-45)  and  (2-51)  as 

A(n+1)(z)  =  A(n) ( z)  +  gCn+1)B(n+1) (z) 
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(2-53) 


where  B^n+^')(z)  is  the  result  of  the  autoregressive  analysis 
of  the  input  u(k) .  If  both  sides  of  (2-53)  are  multiplied 
by  the  input  U(z)  and  then  transformed  to  the  time  domain  the 
result  is  the  estimate  of  the  output  signal  y(k)  for 

y Cn+D (k)  _  y(n>(k)  +  g(n+l)-(n+l)(k)  (2-54) 


where  e^n+’'"^(k)  is  the  backward  error  signal  of  the  autore¬ 
gressive  analysis  of  u(k) .  The  desired  lattice  structure 
is  shown  in  Figure  2.7. 


Figure  2.7  A  Second  Order  MA  Lattice  Structure 
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Recalling  that  the  backward  error  signals  of  different 


stages  of  an  AR  lattice  are  orthogonal,  it  is  seen  that  the 
MA  estimate  of  y^n^(k)  is  a  linear  weighted  sum  of  these 
orthogonal  signals. 

C.  MULTICHANNEL  MODELING 

The  models  previously  considered  have  been  of  single 
input, single  output  systems.  Multiple  input,  multiple 
output  systems  can  be  modeled  as  vector  generalizations  of 
the  models  heretofore  developed. 

1 .  A  Multichannel  Generalization 

A  single  channel  output  of  a  Q-channel  autoregressive 
model  is  the  weighted  summation  of  the  past  outputs  of  the  Q 
channels  of  the  system  to  be  modeled.  Similarly,  a  single 
channel  output  of  a  Q-channel  moving  average  model  of  a 
system  is  the  weighted  summation  of  the  past  N  inputs _ of  .the 
Q  channels  to  the  system  to  be  modeled.  Either  of  these 
statements  may  be  written  in  equation  form 

^  Q  N 

x . (k)  =  Z  £  d. . (n)  x. (k-n)  (2-55) 

J  i  =  l  n=l  ^  1 


Q  =  number  of  channels 
n  =  time  delay  index 
i  =  input  channel 


k  =  time  index 

j  =  output  channel  being 
estimated 
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g 


Xj(k)  =  estimate  of  output  channel  j  at  time  k 

x^(k-n)  =  input  to  channel  i  at  time  (k-n) 

d..(n)  =  weighting  factor  for  estimating  the  output  signal 
1-'  at  channel  j  at  time  k  due  to  input  x  (k-n) 

As  an  example,  for  a  model  with  3  inputs  (Q  =  3), 

(2-55)  would  be  expanded  as 


N 


N 


N 


x^(k)  =  Z  d^(n)  x^(k-n)+  Z  d2^(n)  x2(k-n)+  Z  d3^(n)  x3(k-n) 
n=l  n=l  n=l 


N 


N 


N 


x2(k)  =  Z  d  2(n)  x^(k-n)+  I  d22(n)  x2(k-n)+  Z  d32(n)  x3(k-n) 
n=l  n=l  n=l 


N 


N 


N 


x3(k)  =  Z  d13(n)  x1(k-n)+  Z  d23(n)  x2(k-n)+  Z  d33(n)  x3(k-n) 
n=l  n=l  n=l 


Equation  (2-55)  can  also  be  written  in  matrix  form. 


x  (k)T=  DT  x  (k)  =  CX  (k)T  D]T 


(2-56) 


where 


x(k) 


T 

X  (k) 


[x1(k)  x2(k)  ...  x^(k)]  (Qxl) 
txj(k)  x2(k)  ...  xj(k)]  (NQxl) 


xT(k) 


Cx^(k-l)  x^(k-2)  ...  x^(k-N)]  (Nxl) 
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D 


dll 

-12 

"  *  dlN 

—21 

• 

-22 

' *  *  d2N 

L 

2 

t 

1 

' ‘ *  dQN 

(NxNQ) 


Cdij(D 


di3(2) 


d. .  (N)] 

13 


(Nxl) 


The  error  of  such  an  estimate  is  the  vector  difference 
in  the  actual  output  vector  and  the  predicted  output  vector. 

e(k)  =  x ( k )  -  x(k)  =  [I  -DT]  X(k) 

The  prediction  error  covariance  matrix  is  a  (QxQ)  matrix 

P  =  e[e(k)  eT(k)]  (2-57) 


If  the  trace  of  the  prediction  error  covariance  matrix  is 
minimized  the  result  is 


RD  =  r 


(2-58) 


where 


T 

ft  =  e  {X  (k)  X  (k)} 


=  e 


>c  i  ( k  )" 


x^(k) 


[^(k) 


. .  x^(k)  ] 


Rx  x 

-l-i 


■  R  T 

— 1— Q 


R  T  R  T 

_ v  x  •  •  •  — v  y 

— Q— 1  ^-Q^Q 


(2-59) 


is  a  matrix  of  autocorrelation  and  cross  correlation  matrices 
and 


x,  (k) 

r  ..  . .  .  r 

r  =  e 

11 X 

• 

• 

• 

[x^(k)  ...  x^(k) ]  = 

-1X1  -1XQ 

•  • 

•  • 

•  • 

x«(k) 

r  .  .  .  r 

Lr3  j 

_Yl  ^  Q_ 

(2-60) 

is  a  cross  correlation  matrix  of  sict  (!JQxN)  and  elements 


rx  x  =  £ 

xi(k-l) 

x^(k)  = 

R 

x.x, 
*  j 

(1) 

-ixj 

x. (k-N) 

(N) 

1  X 

Li  5 

— 

(2-61) 
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As  imposing  as  this  may  seem,  it  is  no  more  than  a 
vector  generalization  of  (2-12)  from  the  AR  problem  and 
(2-40)  from  the  HA  problem.  Each  correlation  coefficient  is 
replaced  with  a  correlation  matrix  in  the  multichannel 
generalizat ion . 

It  will  be  useful  in  the  sequel  if  the  two  channel 
lattice  is  examined  additionally.  Equation  (2-56)  expanded 
for  two  channels  is 


x1(k) 

T 

111 

*21 

x-^(k) 

=  D  X  (k)  = 

x2  (k) 

*12 

— 22_ 

x2  (k) 

N 

dn(n) 

d21(n) 

x^(k-n) 

Z 

n=l 

d12(n) 

d22 (n) 

x2 (k-n) 

The  two  channel  error  is  the  vector  difference  in  the  signal 
to  be  estimated  and  the  estimation. 


e  (k)  = 


x^(k) 

x^(k) 

(x^(k)-x^(k) ) 

x2  (k) 

x2(k) 

(x2(k)-x2(k) ) 

(2-63) 
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Equation  (2-58),  the  result  of  a  minimization  of  the  error 
equations,  is  in  two  channels. 


V x„ ( -1) 


Important  to  the  solution  of  (2-58)  is  the  form  of 

the  (NQxNQ)  correlation  matrix  R.  The  matrices  on  the  main 

T 

diagonal  are  all  correlation  matrices  of  the  form  R  ant 

i— i 

so  are  symmetric  and  Toeplitz.  The  off  diagonal 
matrices  are  cross  correlation  matrices  that  are  Toeplitz 
but  not  symmetric.  Symmetrically  opposed  matrices  are, 
however,  transposes  of  one  another,  with  the  result  that  the 
matrix  R  is  also  symmetric.  This  block  Toeplitz  form  is 
amenable  to  a  block  form  of  Levinson  algorithm  solution. 

The  vector  generalization  of  the  AR  lattice  is  shown 
in  Figure  2.8  and  is  described  by  the  equations 

e(n+1)(k)  =  e(n)(k)  -  K(n+1)  e(n)(k-l)  (2-65) 

e(n+1)(k)  =  e(n)(k-l)  -  K(n+1)  e(n)(k)  .  (2-66) 

and  the  initial  conditions 

e(0)(k)  =  e(0)(k)  =  x(k)  (2-67) 

The  signal  paths  of  the  single  channel  AR  model  have  been 
replaced  with  parallel  vector  signal  paths,  adders  with 
vector  adders  and  multipliers  with  matrix  multipliers. 


Figure  2.8  A  Two  Stage  Multichannel  AR  Lattice 

Lattice  predictor  coefficients  K  and  K  may  be 
also  found  from  vector  generalizations  of  the  single  channel 
case  using 

-1  T 

K(n+1)  _  p-(n)  A<n)  (2-68) 


As  in  the  single  channel  case,  the  optimum  predictor 


for  all  orders  less  than  N  are  found  with  the  N-th  order 
solution.  Also,  orthogonality  applies  to  the  backward 
channel  error  signals 

.  _.T 

e  {eik)  e-*  (k) }  = 

and  successive  stages  are  thereby  decoupled. 

For  a  multichannel  MA  model  equations  (2-65)  to  (2-72) 
are  used  with  a  vector  generalization  of  (2-54). 

y;(n+1)(k)  =  X(n)<k)  +  G(n+1)Te(n+1)(k)  (2-73) 

The  error  of  the  MA  model  is  expressed  by 

e^n\k)  =  ;y(k)  -  y  (k) 

and  can  be  minimized  with  the  result 

Q(n)  .  p(n)-1  e{e(n)(k)  vT(k)}  (2-74) 

Figure  2.9  illustrates  the  resultant  structure. 

In  the  next  chapter  the  multichannel  generalization 


or 


p(i) 

V  L 


i  =  3 


is  applied  to  ARMA  modeling. 


III.  THE  AUTOREGRESSIVE  MOVING  AVERAGE  MODEL 


IN  TIME  AND  FREQUENCY  DOMAIN 

While  the  MA  and  AR  models  are  suitable  for  many  appli¬ 
cations,  modeling  some  systems  can  require  a  very  high 
order  model  to  obtain  suitable  results-  Consider  the 
geometric  series  expansion. 

oo 

1  --  n  — n  I  —1 1  , 

- — -I —  =  £  a  z  |az  |  <  1 

(1  -  az”1)  n=0 

This  shows  that  a  single  pole  is  equal  to  an  infinite  number 
of  zeros,  and  conversly  a  single  zero  is  equal  to  an 
infinite  number  of  poles.  Thus  a  system  with  a  single  pole 
will  require  an  AR  model  of  infinite  order,  as  would  a 
system  with  a  single  zero  require  an  infinite  order  MA 
model,  to  be  precisely  represented.  The  ARMA  model  emulates 
both  poles  and  zeros,  and  so  can  model  pole-zero  systems 
without  resorting  to  infinite  order  solutions. 

This  chapter  reviews  contemporary  ARMA  lattice  struc¬ 
tures,  then  develops  some  new  formulations  using  frequency 
domain  data  and  considers  their  implementation  as  a 
performance  evaluation  tool. 
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A.  THE  ARMA  EQUATIONS 

To  model  a  system  with  both  poles  and  zeros  it  is 
intuitively  appealing  to  construct  a  model  that  also  has 
both  poles  and  zeros.  As  was  earlier  discussed,  the  equation 
error  model  will  result  in  linear  equations  as  the  MMSE 
solution.  If  the  condition  that  the  order  of  poles  equal 
zeros  is  imposed  the  equation  error  model  output  estimate  is 
described  by 

.  N  N 

v(k)  =  I  a  u(k-n)  +  Z  b  y(k-n)  (3-1) 

n=0  n  n=l  n 

T  I  t  — 

=  (k) | u1 (k) ]  - 

a 

where  vectors  ^(k),  u(k) ,  b  and  a  are  defined  and  dimensioned 
as  in  the  AR  and  MA  discussions  of  Sections  II. A  and  II. 3 
and  y(k)  is  the  ARMA  output  estimate.  This  constraint  is 
not  limiting  since  some  of  the  coefficients  of  the  numerator 
or  demoninator  can  be  zero.  Minimizing  the  mean  squared 
error  with  respect  to  the  parameter  vectors  a  and  b  yields 
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with  a  MMSE  of 


En,in  =  Syy<0) 


^Tcpt!iToPt] 
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(  3-3) 


Equation  (3-2)  can  be  readily  solved  for  b  ^  and  a  ^  by 

J  —opt  -opt  J 

matrix  manipulation  but  again  a  recursive  in  order  solution 
can  be  derived. 

3.  A  RECURSIVE  IN  ORDER  ARMA  SOLUTION 

The  ( ( 2N+l)x( 2N+1) )  correlation  matrix  of  (3-2)  is 


neither  Toeplitz  nor  symmetric,  but  is  nearly  so.  The 
offending  elements  are  those  associated  with  the  a.  term. 


If  it  is  assumed  that  a  can  be  found  by  other  means,  (3-2) 

o 

can  be  restructured  so  that  the  correlation  matrix  is  block 
Toeplitz  by  subtracting  the  aQ  terms  from  both  sides  of  the 
matrix  equations  (3-2). 
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(3-4) 


Since  aQ  is  just  the  DC  gain,  and  can  be  easily  calculated 
this  is  a  reasonable  assumption.  Equation  (3-4)  can  be 
rewritten  using  the  notation  of  (2-64),  as  was  used  in  the 
discussion  of  the  two  channel  AR  lattice  in  Section  II.C 
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-yy  -yu 
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(3-5) 


where  the  vector  a  contains  the  N  elements  a^...aN  of  the 

T 

vector  a.  Here  use  has  been  made  of  the  identity  R  =R  . 
—  J  — u  v  — vu 

This  is  compared  with  (2-64)  with  x1  =  y  and  x„  =  u 
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To  solve  (3-5)  using  the  algorithms  for  the  two  channel 
lattice  it  is  only  necessary  to  make  the  substitution 


(3-6) 


after  solving  for  the  vectors  d^.  in  order  to  find  b  and  a. 
Since  a^  is  the  DC  gain  of  the  system  to  be  modeled  it  can  be 
found  using 


_  s  (y(k)u(k) } 
0  ”  e  (u(k)u(kH 


(3-7) 


1 .  An  ARMA  Lattice  Structure 

Since  the  equation  error  ARMA  model  can  be  solved 
in  terms  of  a  two  channel  AR  recursive  formulation,  it  is 
not  difficult  to  modify  the  AR  two  channel  lattice  to 
represent  the  ARMA  system. 

The  ARMA  error  can  be  written  as  the  difference 

A 

between  the  output  y(k)  and  the  estimated  output  y(k)  or 
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eQ(k)  =  y(k)-[^‘(k) | uT(k)] 


(3-8) 


from  the  error  equation  model  of  Figure  1.3  and  (3-1).  The 
AR  error  can  be  written  by  substituting  (2-62)  into  (2-63), 
with  y(k)  for  x^(k)  and  u(k)  for  X2(k),  and  transposing. 


so  that 


e(k) 


«T* 

V  (k)l  rp  |  T  I  “t  l 

=  u(k)j  -  ^(k)Tsu«(k)T]  ‘-11  - 


-11  ^12] 

-21  —2  2j 


(3-9) 


=  [e^(k)  e^(k) ] 


where  u'(k)  is  an  N  element  vector  containing  the  elements 
u(k-l) . . . u(k-n)  of  the  vector  u(k). 

If  (3-9)  is  post  multiplied  by  a  vector  tp  where 


-U 


(3-10) 


the  result  is 


e  (k)<2/  =  C e  (k)  e  (k)  ] 
—  y  u 


[:,] 


y(k)-aQu(k)  -  [^(k)  u'T(k)]  £---j  = 


e  (k) 
0 


(3-11) 
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which  is  exactly  the  ARMA  error  of  (3-8),  Furthermore,  the 
mean  squared  ARMA  error  is 


ete^Ck)}  =  ^  e{e(k)  e(k)^|-  ^ 
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(3-12) 


where  P  is  the  forward  prediction  error  covariance  matrix  of 
the  two  channel  AR  model.  If  (3-12)  is  minimized  with 
respect  to  a&  the  result  is 


e  (eu(k)ey(k)} 
e  {eu(k)l 


(3-13) 


which  is  an  alternative  solution  for  aQ. 

To  construct  an  ARMA  lattice  it  is  only  necessary  to 
derive  the  ARMA  error  using  (3-11)  from  the  AR  lattice. 

A  second  order  lattice  model  is  shown  as  Figure  3.1, 

From  (2-65),  forward  prediction  in  the  two  channel  AR 
lattice  can  be  written 

(n) 
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*12 

k2  2_ 

60 


— (  2  ) 


From  this  equation  e^  ;(k)  and  e  (k)  can  be  obtained 


e(n)(k> 


=  e 


(nn)Ck)  ♦  k<n+1>e(n)(k-l)+k'j+1,;'n,<k-l) 


'll  y 


(3-15) 


e(n) (k) 
u 


e(n+l)(k)  +  k55+1)e(n)  (k-l)+k^+1)e(n)  (k-1) 

u  1 2  y  22  u 


(3-16) 


Equations  (3-15)  and  (3-16)  describe  the  synthesis 
structure  that  can  be  used  to  implement  the  ARMA  model  as 
shown  by  Figure  3.2.  For  an  N-th  order  realization  an  input 
is  needed  for  e^^(k).  From  (3-11) 

e*N)  =  eQ(k)  +  aQe<N)(k)  (3-17) 


and  if  the  ARMA  model  is  an  accurate  representation  eQ(k) 
will  be  small  so  that  (3-17)  can  be  rewritten 


(N) 

e 

y 


(k) 


(N)  ... 
a0eu  Ck) 


(3-18) 


This  provides  the  input  to  e^(k)  shown  in  Figure  3.2. 

C.  LATTICE  ALGORITHMS  IN  THE  FREQUENCY  DOMAIN 

In  the  foregoing  it  has  been  tacitly  assumed  that  the 
data  available  for  analysis  is  in  the  form  of  a  finite 
length  time  domain  signal.  This  may  not  be  the  case.  The 
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data  could  be  in  the  form  of  the  Fourier  transform  of  the 
time  domain  signal,  or  it  may  prove  advantageous  to 
process  the  data  in  that  form.  Makhoul  [9]  has  shown  that 
linear  prediction  can  be  done  equivalently  in  either  time 
and  frequency  domain,  but  there  has  been  in  the  literature 
no  exploitation  of  this  for  the  lattice  algorithms.  Thus 
it  seems  useful  to  examine  the  lattice  structures  when  the 
data  is  represented  as  a  frequency  domain  spectrum. 

1 .  The  Single  Channel  A R  Lattice 

The  time  domain  lattice  was  formulated  using  the 
time  domain  equations  (2-24)  and  (2-25).  With  frequency 
domain  data  available  the  single  channel  AR  lattice  can  be 
formulated  using  (2-23)  and  (2-23a).  Its  structure  is  shown 
as  Figure  3.3.  A  unit  time  delay  becomes  a  multiplication 
by  e'3w  with  the  input  signal  represented  by  the  continuous 
spectrum  Y(uj).  Note  that  in  the  backward  channel  the  signal 
E^n^(ai)  is  a  delayed  version  of  E’^^w).  The  frequency 
domain  error  of  a  zero  order  prediction  is 


.(0) 


(u)  =  E’(0)(u>)  =  Y(w) 


(3-19) 


From  the  forward  error  channel  recursion  (2-23)  the  magnitude 
squared  error  of  order  (n+1)  can  be  written  as 

it  it 

-(n+1),  .  r(n+l)/  \  u(n+l)=-(n) ,  ,irr(n),  > 

E  (w)  E  (w)  =  Li  (ui)-k  E  (w)JLE  (o>) 


-k(n+1)E(n)  (U)] 


(3-20) 


Figure  3.3  A  Second  Order  AR  Lattice  in  the  Frequency 

where  E*  indicates  complex  conjugate.  The  total  magnitude 
squared  error  can  be  found  as  the  integral  of  C3-20) 

f  |E(n+1)(w)  |daj  =  /  | E(n)  (w)  |  2d 

—  00  — oo 


■k'nn)t  / 


[E<n)U)f<n) 


(oj)]du/+/  [E 


( n) 


(a>)E^n\w)]dui] 
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k(n.u2 
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—oo 


|E(n)  (oj)  |  2dw 


(3-21) 


Since  XZ*  «  [X*Zl*  the  bracketed  term  of  the  right  side  of 
(3-21)  can  be  rewritten  as 


65 


Lf 


[E(n)(w)E<n) 


gO  « 

(w)]dw  +  /  CE(n)  Ca))E(n)(w)]dw] 


*  * 


=  /  [E(n)(a))E(n)  (u>)]dw  +  /  [E(n)  (ui)E(n)  («)]  dw 

GO  00 

(3-22) 


Since  Y(ui)  is  the  spectrum  of  a  stable  and  causal  sequence 
then  Re{Y(u»)}  is  an  even  function  and  Im{Y(u>)}  is  an  odd 
function.  This  is  also  true  for  the  error  signals  E^n^(u) 
and  E^n^(u>),  The  result  of  the  integrations  of  (3-22),  must 
therefore  be  real  and  the  two  terms  of  (3-22)  are  equal,  or 


/  [E(n)(w)E(n)  (w)Ildu  +  / 


[E(n)  (w)E(n) 


(to)  ]dw 


oo  *»  oo  A 

=  2/  [E<n)(ui)E(n)  (u)]dw  =  2/  [E(n)  (aj)E(n)(u))]du 


(3-23) 


Substituting  (3-23)  for  the  middle  term  of  (3-21) 

(n+1 ) 


and  minimizing  the  result  with  respect  to  k 

sV 


yields 
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/  CE(n)(uj)E(n)  (uj)  ]da> 


/  |E(n) (u) |2dw 


for  a  continuous 
spectrum  Y(u>)  (3-24) 


—  CO 


66 


or  equivalently 


.(n+1)  _  —  tt 


f  [E^n\w)E^n^  Cu))]doj  for  a  periodic 


(n) ,  >  1 2, 


/  iriw(U) rdu 


spectrum  Y(w)  (3-25) 

with  period  2n 
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.(n+1)  k=0  mT  mT 


V  \t"hy\2 

k=0 


for  a  discrete 
periodic  m-point 
spectrum  v,k  .  (3-26) 

1  '•mrr 

with  period  mT 


If  the  discrete  correlation  theorem  and  Parseval’s 
theorem  are  applied  to  (3-26)  the  result  is 


,.(n+l)  _  e(e(n)  (k)e(n)(k-l)}  _  (n+1) 

K  -  - -J- -  -  <F 

e{e(n)  (k)} 


( 3-27) 


which  is  equivalent  to  the  forward  solution  (2-27). 

The  same  procedure  can  be  applied  to  the  backward 
error  recursion  (2-23a)  with  the  result 


00 

/  [E(n)(w)E(n)  (u))]dw 
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for  a  continuous 
spectrum  Y(u>)  (3-28) 
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for  a  periodic 
spectrum  Y(w)  (3-29) 
with  period  2ir 
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(3-30) 


which  can  be  shown  to  be  equivalent  to  (2-28),  the  backward 
equation . 


k(n+D 


e {e(n)(k)e(n) (k-1)} 
e(e(n)  (k)} 


.(n+1) 


(3-31) 


2 .  The  Two  Channel  ARMA  Lattice 

The  multichannel  lattice  can  be  similarly  developed 
using  frequency  domain  transformations  of  the  time  domain 
equations.  Doing  so,  the  two  channel  lattice  equations 
become 


i 
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(3-33) 
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T 

=-,(n+l),  _  =(n),  7(n+l)  -(n),  .  ,,  „  n  % 

E'  (u>)  =  E  (ui)-K  E  C  a) )  (3-40) 

where  the  limits  of  integration  depend  upon  the  periodicity 
of  the  available  data.  If  the  data  is  in  a  discrete  form 
the  integrations  are  replaced  with  summations. 

It  is  important  to  note  that  when  the  data  is  in  a 
discrete  form  the  result  of  a  correlation  obtained  for  a 
time  shift  other  than  zero  will  differ  when  calculated  in 
the  frequency  domain  from  the  result  normally  obtained  from 
time  domain  data.  Consider  first  two  data  signals  y(k)  and 
x(k),  both  defined  for  l_<k_<N,  and  k  an  integer.  The 
correlation  of  delay  l  is  obtained  by  the  summation 

"N-?. 

R„(i)  =  £  y(k)  x  (k+£ 

yx  _k=l 
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and  has  the  effect  of  windowing  the  data  so  that  data  points 
outside  the  window  are  ignored.  If  these  data  signals  are 
transformed  to  the  discrete  frequency  domain  with  a  digital 
Fourier  transform  the  resultant  spectrum  is  periodic,  which 
implies  that  in  the  time  domain  the  data  signal  has  also 
been  made  to  be  periodic.  A  correlation  obtained  in  the 
frequency  domain  is  now 

R  U)  =  i-/*  YU)  X(-w)  e%^  dw 
yx  2  it  N 

J  —IT 

which  when  transformed  back  to  the  time  domain  is  equivalent 
to 

N-i  N  1  , 

R  Al).  =  Z  y (k)x(k+£ )  +  Z  y(k)x(k+£-N)  rr~ 
yx  k=l  k=N-a  +  l  w"x’ 

This  is  because  the  effect  in  the  time  domain  of  the  fre¬ 
quency  domain  time  delay  eij-jp  is  a  circular  shift  of  the 
data  in  the  time  domain.  As  was  discussed  in  the  pre¬ 
ceding  chapter  time  domain  data  would  normally  use  a  linear 
shift  and  then  be  either  windowed  or  would  be  considered  of 
infinite  length.  Either  method  introduces  it's  own  dis¬ 
tortion  into  the  estimation  of  the  correlation.  The  circular 
shift  of  the  frequency  domain  correlation  is  a  third  and 
different  kind  of  distortion.  This  is  important  when  the 
number  of  data  points  is  small,  because  a  bias  will  be 
introduced  into  the  estimate  of  the  correlation. 


D.  SYSTEM  IDENTIFICATION  WITH  THE  ARMA  LATTICE 

The  problem  of  system  identification  concerns  the 
determination  of  a  mathematical  model  of  a  system  from  its 
input  output  relationship.  An  impetus  for  this  research  has 
been  the  issue  of  fault  detection  within  a  system.  If  a 
model  can  be  constructed  of  an  operating  system,  either  on 
line  or  under  test,  the  model  parameters  might  reveal  if  the 
system  is  functioning  normally  and  if  not,  what  class  of 
component  is  faulty.  While  the  ARMA  lattice  has  charac¬ 
teristics  which  make  it  particularly  attractive  for  this 
application,  particularly  its  robustness,  there  are  some 
characteristics  that  still  deserve  attention. 

In  system  identification  the  lattice  model  would  be 
applied  as  shown  in  Figure  3.4.  The  lattice  parameters  K^n\ 
and  ag  that  result  from  a  test  can  be  compared  with  a 
dictionary  of  parameters  obtained  under  normal  and  faulty 
conditions.  Using  lattice  parameters  directly,  though, 
results  in  more  model  parameters  than  if  the  vectors  a 
and  b  were  used  for  identification.  Identification  using  a 
and  b  uses  (2n+I)  parameters  for  an  n-th  order  system  while 
lattice  parameters  result  in  (8n+l)  parameters.  Therefore 
it  is  beneficial  to  reduce  the  number  of  lattice  parameters 
that  must  be  examined  in  the  identification  problem. 
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Figure  3.4  ARMA  Lattice  Used  for  System  Identification 

1 .  Restraints  on  ARMA  Lattice  Parameters 

ARMA  lattice  parameters  are  a  function  of  input 
signal  second  order  statistics.  The  input  signals  u(k)  and 
y(k)  to  the  ARMA  lattice  are  related  by  the  transfer  function 
of  the  system  under  test.  This  suggests  that  the  ARMA 
lattice  parameters  are  interrelated.  To  examine  the  pro¬ 
blem  the  equations  for  the  lattice  parameters  will  be 
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expanded  in  terms  of  the  system  input  u(k)  and  output  y(k). 
The  vector  defining  the  initial  conditions  for  the  ARMA 
lattice  first  stage  is 


(3-41) 


The  forward  and  backward  prediction  error  covariance 
matrix  become 


s[e(0)(k)e(0) 


(k)] 


y2(k) 
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(  0 ) 

and  the  matrix  A  becomes 


A(0)  =  e[e(0)(k)e(0)(k-l)] 
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y(k)y(k-l)  y(k)u(k-l) 
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(3-43) 


Now  consider  the  result  when  the  signals  u(k)  and 


y(k)  are  from  a  system  as  shown  in  Figure  3.5 
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Figure  3.5  System  Excited  with  White  Noise 


From  linear  system  analysis  theory  several  statements 
can  be  made  about  the  second  order  statistics  of  y(k)  and 
u(k),  assuming  that  u(k)  is  white  noise. 
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These  identities  can  be  substituted  into  (3-42)  and 
( 3-43 ) •  Thus 
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With  these,  the  lattice  parameters  can  be  found 
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From  (3-51)  and  (3-52)  some  restraints  on  the 
lattice  parameters  are  evident. 
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It  can  thus  be  seen  that  it  is  not  necessary  to 
calculate  all  eight  lattice  parameters.  Only  three,  plus 
the  DC  gain  aQ ,  are  needed.  Alternatively,  the  parameter 
matrices  can  be  calculated  and  either  (3-54)  or  (3-55)  used 
for  an  alternative  equation  for  calculation  of  aQ . 

Equation  (3-53)  can  be  interpreted  in  light  of  the 
lattice  structure,  k^.^  and  k^^  are  the  coefficients  that 
are  the  predictors  of  the  future  value  of  u(k) .  Since  u(k) 
is  uncorrelated  from  sample  to  sample  no  prediction  can  be 
made  and  these  must  be  zero. 

Some  of  these  relations  can  be  extended  to  later 
stages.  The  forward  error  vector  of  the  first  stage  is 
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because  of  (3<-53)  while  the  backward  error  vector  is 
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where  F^{*}  indicates  a  linear  function  of  the  arguments. 
A^^  is  novj  defined  as 
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because  of  the  independence  of  samples  of  u(k) .  When  this 

(  2 ) 

is  applied  to  (2-65)  for  k  the  result  is 
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By  induction  this  argument  can  be  extended  to 


arbitrary  order  with  the  result 
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Using  (3-61)  the  forward  error  vector  can  be  written 
for  order  n 
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and  this  can  be  used  to  find  the  forward  error  covariance 
matrix 
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If  this  is  now  used  with  (3-59)  of  order  n  to  find 
the  backward  Dararaeter  matrix  the  result  is 
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and  from  this  it  is  seen  that 
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For  the  single  channel  AR  lattice  the  forward  and 
backward  coefficients  were  equal,  but  for  the  multichannel 
lattice  it  is  known  that  the  forward  and  backward  parameter 
matrices  are  not  equal.  Nuttal  [7]  has  shown  that  the 
determinants  of  the  forward  and  backward  error  covariance 
matrices  are  equal,  and  from  this  it  can  be  shown  that  the 
determinants  of  forward  and  backward  lattice  parameter 
matrices  are  equal. 
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It  is  known  that  P 


so  it  follows  that 


(0)  _  p(0) 

Det(P(0)) =  Det(P(0))  (3-67) 

Also 

p(n+l)  _  p(n)  ^  n+1)  K(n+l)-j 

_  p(n)  g(n+l)  ^(n+1)  1  _  ^(n+l)-j  (3-68) 

and 

p(n+l)  _  p(n)  j-j  ^  K(n+1)  ^(n+l)j 

_  p(n)  ^(n+1)  _^(n+l)-j  ^(n+1)  (3-69) 

It  is  known  that 

Det(AB)  -  Det(A)Det(3)  =  Det(3)Det(A) 
so  that  the  determinants  of  (3-68)  and  (3-69)  can  be  written 

Det(P(n+1))  =  Det(P(n)  )Det(K(n  +  1))Det(K(n+1)  -K(n+1)) 

Det(P(n+1))  =  Det(F(n)  )Det(K(n+1))Det  (K(n+1)  -K(n+1)) 
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which  shows  that  Det(P(n+1))  =  Det(F(n+1))  if  Det(P(n) )=Det(F(n) ) . 
In  view  of  (3-67)  it  can  be  inductively  reasoned  that 

Det(P(n))  =  Det(FCn))  (3-70) 

Now  the  determinats  of  the  lattice  parameters  defined  by 
(2-63)  and  (2-69)  can  be  written 

Det(K(n+1) )  =  Det(P(n)  1)Det(A(n)T) 

Det(K(n+1))  =  Det(P(n)  1)Det(A(n)) 

and  since  Det(A)  =  Det(AT) 

Det(K(n+1))  =  Det(K(n+1))  (3-71) 


A  summary  of  these  restraints  is  presented  as  Table  I. 


TABLE  I 

SUMMARY  OF  RESTRAINTS  ON  ARMA  LATTICE  PARAMETERS 


Restraint 


Det  (K 

(n) 

k12 

5-(n) 

,<2i 


T7T7 


=  k 


(n) 

22 

(1) 

11 


)  =  Det(K 
(n) 


L?(n) 


=  0 

rr(n) 
*11 

-a  k(n) 
0  12 


22 

"ao 


-  K 


(1) 

11 


k(1) 
<2  2 


) 


none 


RuuU) 

Ruu(*> 

Ruu(*> 


Ruu(*> 


Conditions 


au5  (JO 
auS  U> 
au5  (JO 
au5  (Jl) 
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2 .  Uniqueness  of  the  Lattice  Parameters 

A  particularly  troublesome  consideration  in  using 
ARMA  lattice  parameters  for  system  identification  is  the 
lack  of  uniqueness  of  the  lattice  parameters  for  the  system 
under  test.  This  is  evident  because  the  lattice  parameters 
are  not  just  functions  of  the  system  under  examination ,  but 
also  of  the  test  signal  used  to  drive  the  system.  Thus  for 
different  test  signals  such  as  an  ensemble  of  finite  white 
noise  input  signals,  an  ensemble  of  lattice  parameters  are 
obtained.  Interestingly,  experiments  have  shown  that  the 
transfer  function  coefficients  calculated  using  these  lattice 
parameters  will  generally  have  less  dispersion  than  the 
lattice  parameters  from  which  they  are  calculated. 

Thus,  the  selection  of  a  test  signal  influences  both  the 
accuracy  of  the  model  and  the  expectation  of  the  lattice 
parameters  that  result.  The  mean  square  value  of  equation 
error  of  the  ARMA  model  can  be  expressed  in  integral  form  as 


E 


1 

7W 


TT 

f 


-TT 


E(e^wT)  E*(e^uT)dwT 


( 3-72a) 


The  equation  error  can  be  written  in  the  z  domain  as 


ECz)  =  [BCz)H(z)  -  A(z)]U(z) 


( 3-72b) 
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Substituting  (3-72b)  into  (3-72a)  results  in 


E  =  /ir|[BCej“T)HCe^T)-A(e;5wT)]!2(LKe;5‘i)T)|2du)T 

■'71 

(3.-73) 

Equation  (3-73)  shows  that  the  power  spectrum  of  the  driving 
source  is  a  weighting  function  for  the  resultant  modeling 
error.  To  obtain  a  model  that  is  equally  accurate  throughout 
the  spectrum,  a  test  source  must  be  used  which  has  equal 
energy  throughout  the  spectrum.  An  obvious  source  with  a 
uniform  power  spectrum  is  a  white  noise  source.  Besides 
giving  equal  weighting  to  the  error  spectrum,  white  noise 
involves,  as  has  been  shown  in  the  preceeding  section, 
simplified  calculations  and  a  simplified  lattice  structure. 
This,  however,  is  not  the  total  answer.  When  applied  to 
system  identification,  the  correlations  and  cross  correlations 
of  y(k)  and  u(k)  generally  will  not  be  available.  Instead 
these  quantities  will  be  estimated  from  finite  sample 
functions  of  y(k)  and  u(k) .  Even  though  these  signals  may 
be  ergodic,  the  values  obtained  are  only  estimates  of  the 
correlations . 

Experiments  with  a  white  noise  source  have  shown 
that  different  ensembles  of  input  signals  will  result  in 
different  lattice  parameters  for  each  member  of  the  ensemble. 
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In  order  to  improve  the  ensemble  averages  obtained  for  the 
lattice  parameters,  some  form  of  input  preprocessing  can  be 
performed  to  whiten  the  input  or  equivalently  flatten  it's 
spectrum. 

Consider  first  Figure  3.6a,  which  is  a  modification 
of  the  system  of  Figure  3.4.  An  adjustable  filter  with 
transfer  function  H(u>)  has  been  added  to  the  input  of  the 
system  under  test.  The  transfer  function  H(w)  can  be  adjusted 
to  made  U'CuO  white  if  the  spectrum  of  the  noise  source  is  not 
white  and  has  not  zeros  on  the  u  axis .  If  the  system  under 
test  is  linear,  then  the  lattice  filter  derived  from  Figure 
3.6b  is  equivalent  to  the  one  derived  from  Figure  3.6c. 

Figure  3.7  is  a  block  diagram  for  calculating  the  frequency 
domain  ARMA  lattice  model  with  four  preprocessing  techniques. 

a)  Preprocessing  in  the  time  domain  on  both  the  input 
and  output  signals  of  the  system  under  test. 

b)  Preprocessing  in  the  frequency  domain  on  both  the 
input  and  output  signals  of  the  system  under  test. 
Since  time  frequency  domain  operations  are  equivalent 

they  can  be  interchanged.  For  comparison  they  are  considered 
separately  to  determine  if  there  are  any  computational 
advantages  of  one  over  the  other.  It  should  be  noted  that 
as  long  as  equivalent  operations  are  performed  on  both  input 
and  output,  as  discussed  earlier,  the  resultant  lattice  model 
will  still  model  the  system  under  test,  though  the  specific 
lattice  coefficients  may  differ. 
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(a) 


Figure  3.6  Lattice  Input  Preprocessing 


8 


Figure  3.6  Lattice  Input  Preprocessing  (con't) 


► 


Lattice  Parameters 


Figure  3.7  Frequency  Domain  ARMA  Lattice  with  Preprocessing 
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A  random  noise  generator  can  be  represented  as  shown 
in  Figure  3.6c.  If  this  generator  provides  the  source  signal 
of  Figure  3.7,  then  preprocessing  to  make  U'(w)  white  is 
equivalent  to  the  multiplication  in  the  frequency  domain  of 
U(id)  by  G~^(w)  ,  assuming  that  the  inverse  exists  and  is  stable 
In  the  discrete  n^point  frequency  domain  U(oj)  can  be 
written  as  a  vector 

UT  =  [U0  U1,  Un_13  ( 3-74a) 

G(w)  in  its  discrete  form  can  also  be  written  as  a  vector 

GT  =  [Gn  G,  ...  G-]  (3-74b) 

—  01  n-1 

Multiplying  by  G  ^(oj)  can  be  expressed  in  the  discrete 

frequency  domain  as  a  multiplication  of  the  vector  U  by  the 
.  T  T 

square  matrix  H,  namely  U'  =U  H,  where 

H  =  [h^j]  and  hij=0  f°r  ^  ^  3  (3-75a) 

and 

h..  =  O’”1,  for  i  =  1,  2,  ...,  n  (3-75b) 

ii  i^-l  7 

where  tL  is  the  expected  value  of  the  ith  harmonic  of  the 
input  signal  U(w)  taken  over  an  ensemble  of  the  input  signals. 


For  H  to  be  useful  GU)  should  have  no  poles  or  zeros 
on  the  unit  circle.  GU)  can  be  known  apriori  (or  estimated) 
from  an  autoregressive  analysis  of  the  signal  generator 
output,  or  by  averaging  the  spectra  of  several  (.possible 
overlapping)  segments  of  the  generator  output.  Equations 
(3-7 5a)  and  (3-75b)  are  useful  when  dealing  with  an  ensemble 
of  input  signals.  Alternately  for  a  single  member  of  the 
ensemble 


h^  =  ^or  *  =  2,  •  ‘  ‘  n  (3-75c) 

can  be  used.  This  is  equivalent  in  the  time  domain  to 
exciting  the  system  with  a  unit  impulse. 

From  Figure  3.6c  UU)  =  NU)GU) .  From  Figure  3.6a 
the  inputs  to  the  lattice  filter  are  given  by 

U'U)  =  UU)HU)  =  NU)GU)HU)  =  NU) 

when 

HU)  =  G-1U)  .  Also, 

Y’U)  =  U’(w)TU)  =  NU)TU) 


and  the  transfer  function  derived  by  the  lattice  model  is 


Y*U) 

U'U) 


TU) 
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From  Figure  3.6b  U» <w)=U<u)H(w)  and  Y' (w)=Y(w)H(w)*U(w)T(w)H(w) . 
The  transfer  function  for  the  lattice  is 


Y'  (w) 
U'  (w) 


T(w) 


which  is  equivalent  to  that  of  Figure  3.6a.  In  Figure  3.6b 
the  spectra  presented  to  the  lattice  filter  are  not  functions 
of  G(uj)  ,  and  the  coefficients  obtained  can  be  used  directly 
to  identify  the  system  under  test. 


E.  FREQUENCY  DOMAIN  ARMA  LATTICE  SIMULATIONS 

To  verify  some  of  these  concepts  a  FORTRAN  program  was 
written  that  implements  the  ARMA  frequency  domain  lattice 
with  preprocessing  introduced  as  shown  in  Figure  3.8.  A 
Hamming  window  can  be  selected  for  preprocessing  in  the  time 
domain.  The  selectable  frequency  domain  processing  imple¬ 
mented  is  that  of  (3-75c).  The  spectrum  Y(w)  is  divided  by 
the  spectrum  U(u)  to  give  Yf(u>),  the  new  system  output 
spectrum.  The  input  spectrum  UCw)  is  replaced  with  a  unit 
magnitude  real  spectrum. 

When  preprocessing  is  not  desired  then  the  inputs  to 
the  ARMA  frequency  domain  lattice  are  changed  to  U'(n)= 

U(n)  and  Y* (n)=Y(n) . 

For  comparison  a  program  was  written  that  solved  the 
ARMA  normal  equations  by  Gaussian  elimination.  Listings 
of  these  progsams  are  available  from  the  author. 
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u(n) 


y  (n) 


Figure  3.8  ARMA  Lattice  Simulation  Program  Block 
Diagram  with  Input  Preprocessing 
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More  than  twenty  different  systems,  of  up  to  fourth 
order,  were  used  as  the  system  under  test.  The  number  of 
data  points  used  for  a  single  analysis  varied  from  32  to 
4096.  The  number  of  points  was  always  an  integral  power  of 
two  because  of  an  FFT  routine  used  in  the  frequency  domain 
lattice . 

To  illustrate  the  ability  of  the  frequency  domain 
lattice  to  accurately  model  a  system  two  examples,  typical 
of  all  the  problems  exercised  are  shown  here.  System  1, 


TABLE  II 
SYSTEM  1 


Transfer  Function  Coefficients 


Numerator  Denominator 


A(  0 )  = 

0.25 

3(1)  =  -1.14 

A(  1 )  = 

0.35 

B( 2 )  =  1.4549 

A(  2  )  = 

0.245 

B ( 3  )  =  -0.849 

A(  3  )  = 

0.0 

B(  4  )  =  C. 40745 

A(  4 )  = 

0.0 

Zeros 

Poles 

Real 

Imaginary 

Real  Imaginary 

-0.70 

0.70 

0.50  0.50 

-0.70 

-0.70 

0.50  -0.50 
0.07  0.90 
0.07  -0.90 

characterized  in  Table  II,  is  a  fourth  order  low  pass  system. 
Figure  3.9a  and  3.9b  show  the  magnitude  of  the  transfer 
function  of  the  model  obtained  by  exciting  the  system  with  a 
single  (256  point)  sample  of  white  noise  from  the  Gaussian 
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elements  and  the  frequency  domain  lattice  using  no  pre¬ 
processing.  As  has  been  observed  in  the  past  [15]  the 
lattice  results  in  a  closer  model  of  the  system.  Figure  3.9c 
and  3,9d  show  the  pole  zero  plots  of  the  resultant  transfer 
functions . 

Figure  3.10a  and  3.10b  are  the  transfer  function  magnitude 
found  with  the  same  methods  using  4096  data  points.  Here  the 
lattice  model  is  still  a  more  accurate  model  than  the  batch 
method. 

While  in  most  cases  the  lattice  offered  superior  perfor¬ 
mance  to  the  Gaussian  elimination  method.  System  2,  Table  III, 
is  a  counter  example  that  shows  that  this  is  not  always  the 
case.  In  this  example  the  system  is  a  second  order  band  pass 
filter.  Figure  3.11a  and  3.11b  are  the  transfer  function 
magnitude  plots  found  using  128  points,  and  Figure  3.11c  and 
3. lid  show  the  pole  zero  plots.  Here  it  is  seen  that  the 
batch  method  produces  a  superior  model.  Figure  3.12  shows 
the  same  system  analysis  using  1024  points,  and  both  methods 
accurately  model  the  system. 

To  determine  the  ability  of  the  input  preprocessing  to 
compensate  for  nonwhite  input  to  the  system  the  arrangement 
of  Figure  3.13  was  used  to  generate  a  colored  input  signal. 

The  plant  used  is  that  of  System  1.  The  test  input  is  a 
white  Gaussian  noise  source  followed  by  a  single  pole  filter, 
with  the  pole  position  varied  from  -0.5  to  0.5  for  different 
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Figure  3.10c  Gaussian  Elimination  4096  Points 
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Figure  3.10d  Lattice  Model  4096  Points 
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Figure  3.11c  Gaussian  Elimination  128  Point 


Figure  3. lid  Lattice  Method  128  Points 
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Figure  3.12b  Lattice  Method  1024  Points 
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Figure  3.12c 
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Figure  3.12d  Lattice  Method  1024  Points 
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Figure  3.13  Input  Preprocessing  Test 


TABLE  III 
SYSTEM  2 

Transfer  System  Coefficients 
Numerator  Denominator 


A( 0 )  =  1,00 
A(l)  =  -1.40 
A( 2 )  =  0.95 


BCD  =  0.40 
B(2)  =  0.89 


Zeros 

Real 

0.700 

0.700 


Imaginary  Real 

0,67823  0.200 

-0.67823  0,200 


Poles 

Imaginary 

0.92195 
-Q. 92195 
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test  runs.  The  filter  w^s  operated  under  three  conditions 
of  preprocessing; 

1)  No  preprocessing. 

2)  Inputs  U(w)  and  Y(w)  normalized  as  in  (3-75c) 

3)  Multiplication  of  the  input  data  by  a  Hamming  window 
and  then  normalizing  in  the  frequency  domain  as  in 
(3-7 5c) . 

Figure  3.14  shows  two  typical  plots  of  the  first  stage 
lattice  parameters.  The  horizontal  axis  is  the  pole  position 
b,  of  the  input  filter  of  Figure  3.13.  The  vertical  axis  is 
the  value  of  the  lattice  parameter,  in  this  case  k^^  and 
k22^  for  cases  1)  and  2).  For  all  first  stage  parameters, 
case  3)  results  in  essentially  identical  results  to  case  2). 
It  is  clear  that  the  preprocessing  has  made  the  lattice 
parameters  independent  of  the  input  filter  characteristics. 

Figure  3.15  shows  two  similar  plots  for  second  stage 
parameters  and  ]<12  .  Again  the  preprocessing  is 

effective,  and  case  3)  is  essentially  identical  to  case  2). 

( 3 ) 

Figure  3.16  shows  parameter  k^  .  The  first  plot  is 
that  of  case  1)  and  2).  Here  the  parameter  shows  some 
sensitivity  to  the  pole  position,  though  much  less  than  the 
unnormalized  input  signal.  The  second  plot  is  for  cases 
1)  and  3).  The  windowed  and  compensated  lattice  still  shows 
no  sensitivity  to  the  input  filter  pole  position.  Of  the 
eight  parameters,  this  one  showed  the  most  variation  with  the 
pole  position. 


STAGE  3  a  KF11  STAGE 


Here  the 


Figure  3.17  is  the  fourth  stage  parameter  • 

parameter  calculated  without  windowing  varies  over  an  even 
wider  range  than  the  parameters  calculated  without  any 
preprocessing.  Several  other  of  the  parameters  of  the  fourth 
stage  do  the  same.  The  parameters  of  the  windowed  and  com¬ 
pensated  lattice,  case  3),  are  all  still  insensitive  to  the 
input  pole  position. 

It  is  apparent  from  this  experiment  that  input  signal 
preprocessing  is  effective  in  compensating  for  the  effects 
of  a  nonwhite  input  signal,  and  results  in  a  consistent 
set  of  lattice  parameters  for  identification.  The  pre¬ 
processing  should  consist  of  windowing  of  the  input  data 
(to  eliminate  the  Gibb's  phenomenon)  and  normalization  to 
whiten  the  signal. 

F.  FILTER  SYNTHESIS  WITH  THE  ARMA  FREQUENCY  DOMAIN  LATTICE 

So  far  the  lattice  filters  have  been  developed  as  models 
for  systems  with  measurable  input  and  output.  It  is  possible 
to  define  the  input  signal  to  the  filter  as  a  unit  impulse 
so  that  the  lattice  models  a  desired  response  in  the  frequency 
or  time  domain.  In  this  application  the  recursive  in  order 
nature  of  the  lattice  solution  becomes  of  interest.  The 
lattice  can  be  expanded  in  order  until  a  suitable  fit  of  the 
desired  specification  is  obtained. 
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Consider  what  happens  when  the  frequency  domain  ARMA 
lattice  is  used  to  model  a  system  when  the  excitation  of  that 
system  is  a  unit  impulse.  The  configuration  is  shown  in 
Figure  3.18.  U(n) ,  the  spectrum  of  the  input  to  the  system 
under  test,  is  a  flat  spectrum  of  unit  magnitude,  exhibiting 
the  previously  discussed  desirable  property  of  equal  energy 
across  the  spectrum.  The  output  spectrum  of  the  system  Y(n) 
is  exactly  H(n),  the  discrete  Fourier  transform  of  the 
impulse  response  of  the  system  under  test. 


Figure  3.18  Frequency  Domain  Modeling 


The  lattice  parameters  obtained  as  a  result  of  this 
analysis  could  be  used  to  construct  a  synthesis  model.  This 
synthesis  model  would  exhibit  the  best  fit,  -in  a  least  mean 
square  sense,  that  a  filter  of  its  order  could  obtain.  The 
order  of  the  lattice  can  be  increased  to  obtain  an  arbitrarily 
small  MSE.  When  designing  to  a  specification  it  is  often 
desired  that  some  portions  of  the  frequency  characteristics 
of  the  resultant  filter  be  weighted  more  heavily  than  other. 
This  too  can  be  accomplished  with  the  ARMA  lattice.  U(n)  and 
Y(n)  can  both  be  multiplied  by  a  weighting  function  which 
provides  preemphasis  to  those  portions  of  the  spectrum  for 
which  a  more  accurate  fit  is  desired. 

1 .  Filter  Synthesis  Examples 

To  illustrate  this  capability  it  was  postulated  that 
a  digital  filter  was  desired  with  a  magnitude 

| H(n) |  =  0.495  cos  (3-^  •  n)  +.505,  0<n<1023 

This  desired  magnitude  was  transformed  to  a  minimum  phase 
spectrum  using  a  digital  approximation  of  the  Hilbert  trans¬ 
form  [14].  This  spectrum  was  used  as  the  input  Y(n)  to  the 
frequency  domain  lattice,  with  a  unit  real  spectrum  applied 
as  the  U(n)  input.  The  magnitude  of  the  first  and  second 
order  models  of  this  desired  spectrum  are  shown  in  Figure 
3.19.  The  second  order  model  magnitude  appears  to  match 
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Figure  3.19b  Second  Order  Synthesis  Model 
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Figure  3.19c  Second  Order  Pole-Zero  Plot 
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the  desired  magnitude,  and  as  shown  by  the  pole  zero  plot, 
Figure  3.19c,  is  stable  and  minimum  phase.  The  digital 
transfer  function  of  this  second  order  filter  is 


H(z) 


0 . 302504+0. 495032z_1+0.202512z 
1+0. 00004z_1+0 . 000005 z“2 


As  another  example  the  frequency  domain  lattice  is 
used  to  find  a  digital  equivalent  of  the  analog  low  pass 
circuit  of  Figure  3.20  which  exhibits  40  dB  attenuation  at 
5  kHz.  This  circuit  was  analysed  using  the  SPICE  program, 
a  circuit  analysis  program.  This  program  produced  an  AC 
analysis  of  the  circuit,  showing  the  complex  output  when  the 
circuit  is  driven  by  a  unit  magnitude  signal  for  frequencies 
from  DC  to  10  kHz,  the  equivalent  of  a  Bode  plot  of  the 
circuit  transfer  function. 

This  complex  spectrum  was  applied  to  the  frequency 
domain  lattice  as  in  the  previous  example.  The  third, 
fourth  and  fifth  order  solutions  are  shown  in  Figure  3.21. 

The  transfer  function  of  the  fifth  order  solution  is 

-. 0046 9-0. 026 9 z'1-0 . 310z'2 -0 .03164 z"3-0 . 316z"4+ . 00096 z 
H(z)  =  — - 

l-2.159z-1+2.534z"2-1.675z"3+.6178z"4-0.0979z“5 


Input 


36.9  mH 


48.8  mH 
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ic 
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".0  513  pf  ^ 
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.00444  uf 
"  .119  uf 

Output 
"  .0675  uf 


Figure  3.20  Analog  Low  Pass  Filter 
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The  pole  zero  and  phase  plots  of  the  fifth  order 
solution  are  shown  in  Figure  3.21.  There  is  an  additional 
pole  not  shown  on  the  plot  that  is  located  at  3=-4.5.  .The 
phase  plot.  Figure  3.21e,  shows  that  the  analog  circuit  has 
a  phase  at  the  Nyquist  frequency  of  approximately  it/4.  Of 
course  the  digital  filter  cannot  match  this  phase,  it  must 
be  either  zero  or  ir  at  the  folding  frequency.  It  is 
believed  that  this  discontinuity  in  the  phase  is  the  cause 
of  the  digital  filters  non  minimum  phase  characteristic. 

The  filter  is  stable,  however,  and  does  meet  the  desired 
magnitude  specification. 


Model 


Magnitude 


+  Model  Pole 
o  Model  Zero 
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AUTOREGRESSIVE  LATTICE  PARAMETERS  CALCULATED 


USING  THE  CORRELATION  MATRIX 


In  the  foregoing  chapters  several  methods  of  solution  of 
the  AR  and  ARMA  normal  equations  have  been  discussed.  These 
equations  can  be  solved  by  matrix  methods  but  for  large 
order  problems  this  is  a  computationally  expensive  method. 

If  the  available  data  is  windowed  to  ensure  a  Toeplitz  and 
symmetric  correlation  matrix  the  equations  can  be  more 
efficiently  solved  using  Levinson's  algorithm.  From 
Levinson's  recursion  the  lattice  structure  can  be  imple¬ 
mented  but  requires  more  calculations  because  it  is  necessary 
to  update  prediction  error  signals  at  each  stage.  Mokhoul 
[10]  developed  an  algorithm  that  implements  the  single 
channel  AR  lattice  but  does  not  require  the  generation  of 
these  error  signals,  and  thus  is  more  efficient  than  other 
methods . 

In  this  chapter  a  new  implementation  is  developed  for 
both  the  single  channel  and  the  two  channel  AR  lattice  that 
implements  the  lattice  structure  without  the  expense  of 
calculating  error  signals.  In  the  final  formulation  there 
is  no  requirement  that  the  correlation  matrix  be  Toeplitz 
or  symmetric,  removing  constraints  on  how  data  is  windowed 
and  allowing  solution  of  the  problem  when  the  system  under 
analysis  is  not  stationary.  However,  since  the  formal  proof 


depends  upon  the  lattice  formulation  which  assumes  symmetric 
Toeplitz  matrices,  two  solutions  are  presented,  one 
assuming  Toeplit  symmetric  correlation  matrices ,  the  other 
omitting  this  restriction.  This  algorithm  maintains  the 
lattice  advantages  of  maximum  entropy  and  robustness  that 
have  made  the  lattice  an  attractive  analysis  structure.  It 
requires  less  storage  than  conventional  lattice  implementations 
and  can  be  extended  efficiently- to  multichannel  structures  as 
are  used  for  nonlinear  systems  analysis. 

A.  THE  SINGLE  CHANNEL  AUTOREGRESSIVE  LATTICE 

In  the  single  channel  autoregressive  lattice  structure 
of  Figure  4.1  the  forward  and  backward  error  signals  for 
each  stage  are  updated  using  the  equations  from  Section  II. B. 

e(n+1)(k)  =  e(n) (k)  -  k(n+1)e(n) (k-1)  (4-1) 


e(n+1)(k)  = 


e(n)(k-l)  - 


.( n+1 ) 


(n) ,,  . 
e  (k) 


(4-2) 


The  reflection  factors  k 


(n) 


can  be  found  from  the  backward 


equation  (2-28) 


( n+1) 
B 


e  {e^n)  (k)iU)  (k-l)> 
e  {eCn)(k)e(h)(k)> 


(4-3) 
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Sir.cj  the  quantities  e^n\k)  and  e^n^(k)  are  suras  of  weighted 
present  and  past  values  of  the  input  y(k) ,  the  numerator  and 
denominator  of  (4-3)  can  be  expressed  as  the  weighted  sum  of 
correlations  of  the  input  vector  y.  An  algorithm  based  upon 
this  observation  can  Je  expected  to  be  computationally  more 
efficient  than  previous  algorithms  if  the  correlations  of  the 
input  signal  are  known  since  it  is  not  necessary  to  update 
the  error  sequences  at  each  lattice  stage. 

From  Figure  4.1  the  error  signal  e^n^(k)  can  be  written 
as 


(n) 


(k)  =  [y(k)y(k-D  .  .  .y(k-n)  ] 
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(n)~ 
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(n) 
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(n) 

'n 


(4-4) 


where 
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(n) 


-k(n)  =  i 
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i=l 


i  =  0 


w(n)  =  -k(n) 
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(4-5) 
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or 


e(n)(k)  =  £(n)  (k)w(n)  (4-10) 

where  w^n^  is  the  inverted  version  of  v/n^.  With  this 
notation,  a  delay  of  the  backward  error  sequence  is  repre¬ 
sented  by  a  shift  downward  of  the  elements  of  the  weighting 
vector  w^n\  with  a  zero  shifted  into  the  top,  or 


e(n)(k-l) 


=  Z 


(n+1) 


(k) 


0-11) 


The  weight  vector  associated  with  e^n\k)  or  e^n)(k)  can 
be  found  recursively.  Substituting  (4-6)  and  (4-11),  the 
equations  for  the  signals  e^n\k)  and  e^n^(k-l),  into  the 
forward  error  recursion  (4-1),  yields 


e'^Ck)  =  y<n)Tw(n)  -  k<n*l)  (n»l)T 


0 

(n) 
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_  (n+l)T 

r(n)1 

0 

-  Z 

w 

•  ^  m  •• 

-k(n+1) 

0 

;(n) 

w 

(4-12) 
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Comparing  (4-12)  with  (4-6)  rewritten  for  order  (n+1),  yields 


„<nn) 


k<n+1) 


Similarly, 
backward  error 


substituting  (4-6)  and  (4-11) 
recursion,  results  in 


into  ( 4 


5<n)(k) 


y(n*1)(k) 


k(n*l>2<n) 


(k)w<n) 


Z<n*X)Oc) 
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TnT 
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k(n+D 


Comparing  this  with  (4-11)  gives  the  recursion 
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(n+1) 


k(n+l> 


4-13) 


2 ) ,  the 


4-14) 


4-15) 


.2' 


It  should  be  noted  that  the  intital  conditions  for  the 
weight  vectors  and  are 

w(0)  =  Cl]  w(0)  =  Cl]  (4-16) 


As  an  illustration,  from  Figure  4.1  it  can  be  seen  that 
the  second  stage  weight  vectors  can  be  written  as 
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n 

AO 

_ 1 

w(2)  = 

-k(1)*k(1)K(2) 

-k(2) 

w<2>  = 

-k(1)+k(1)k(2> 

_1 

(4-17) 

Applying  these  to  (4-13)  yields 
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-k(3) 


128 


and  to  (3-15) 


O  1 

_ \ 

i 

** 

»(3)  -- 

-k<2) 

-k(1)+k(1)k(2) 

-k<3) 

-k(1>.ka,k(2> 

_ 
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l _ 

-.k(3) 

-k(2)*ka)k(3)-ka)k<2,k(3) 

-ka>+k'1)k(2,+k<2)k(3) 


(4-20) 


yielding  the  third  stage  weight  vectors. 

The  numerator  of  equation  (4-3)  can  now  be  rewritten 
using  the  definition  of  the  error  signals  in  terms  of  the 
weight  vectors 


e  (e(n)(k)e(n)(k-l)} 
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CS‘“>>  Ryy 
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.0 

TnY 


(4-20) 


It  is  seen  that  the  numerator  of  (4-3),  required  for  calcu¬ 
lation  of  can  be  calculated  with  no  presumption  of 

symmetry  or  requirement  that  the  correlation  matrix  be 
Toeplitz.  The  lattice  method  may  still  be  used  if  the 
signals  under  analysis  are  not  stationary  or  estimates  of 
the  correlation  matrix  are  available  that  do  not  result  in  a 
Toeplitz  and  symmetric  matrix.  However,  if  the  correlation 
matrix  of  (4-20)  is  Toeplitz  and  symmetric  the  matrix 
multiplication  may  be  written  as  summations. 


(e(n)(k)e(n)(k-l)}  •=  R  (0)  E  w.wn  ..1 

yy  i=1  i  n-i+i 


n 


n 


+  R  (1)C  E  w.w  .+  E  w.w  . 
yy  i=0  1  n~  i=2  1  n_1  2 


+  Ryy(2UJ0 
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These  summations  are  recognized  as  convolutions  of  the 
weight  vector  ,  or  the  correlation  of  the  forward  and 
backward  weight  vectors, 

Defining  correlation  functions  of  w^  and  w''  as 
/  )  n+l-i, 

4  <*>  =  t  "i'Vl-t-i’  *=1.2,. ..,nn 
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(4-21) 


$(n)U) 


n 

•ill  wiwn+l-a-i’  ^=n>n+l, .  .  .  ,-l  ,0 

X  —  X  — 


2  Z 
i= 1— A 


w.  w 

i 


n+l-i- 


i. 


£,=n,m-i 


,-1,0 


v  _  nU+1 
n - 2 

2  Z 
i=l -l 


and  ( n+  A )  even 


w.w 

i 


n+l-Jl- 


.  +w 

i 


2 

n-£+l 


7 


%-n ,n+l ,...,-1,0 
and  (n+2.)  odd 


(4-22) 


$(n)(~n-l)  =  0 


C 4-2  3 ) 


Equation  (4-.20)  can  now  be  rewritten  as 
E{eCn)(k)i(n)(k-l)}  = 

Ryy(0)$(n)(0)  +  Ryy(l)C$(n)(l)+$(n)(-l)3 
+  Ryy(2)[$(n) (2)+$(n) (-2)3 

+  ...  + 


+  Ryy(n+l)U(n)(n+l)+<l>(n)(-n-l)3  (4-24) 


This  expression  can  be  written  more  compactly  as  the  product 
of  two  vectors,  one  with  elements  of  Ryy ,  the  other  of  4>. 


e{e(n)(k)e(n)(k-l)}  = 


CRyy(0)'-‘Ryy(n+1)] 


$Cn) (0) 

$(n)(l)+$(n)(-l) 


$(n) Cn+l)+<f(n) ( -n-1 ) 


(4-25) 
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or 


e{e(n)(k)eCn)(k-l)} 


(4-26) 


where 


r(n)  =  [Ryy(0)...Ryy(n*l>] 


and 


$(n) (0) 

<fr(n)(l)+$(n)(-l) 
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$(n)(n+l)+$(n)(-n-l) 


(4-27) 


(4-28) 


The  denominator  of  equation  (4-3)  could  be  found  using 
similar  arguments,  but  generally  offers  no  advantage  to 
previously  described  methods.  The  recursion 

e{e(0)(k)e(0)(k)}  =  Ryy<0)  =  p(0) 

2 

e(e(n)(k)e(n)(k)}  =  P(n*'1)  (l-k(n)  )  =  p(n)  (4-29) 

should  be  used,  Equation  (4-3)  may  now  be  rewritten  using 
(4-20)  or  (4-26)  and  (4-29). 
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k(n+1) 


r(n)T  W(n) 
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For  R  Toeplitz  and 
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(4-30a) 


,  ,T| 

[w(n)  | 0  ]R 

0 

n(n) 

For  general 
R 


(4-30b) 


1 .  Implementation  of  the  Single  Channel  Algorithm 

The  lattice  can  be  implemented  using  correlation 
weights  by  applying  (4-16),  (4-13),  (4-27),  (4-20)  or 
(4-28),  and  (4-29)  and  (4-30).  A  flowchart  of  the  procedure 
is  presented  as  Figure  4.2. 

For  an  mth  order  analysis  the  procedure  will  be: 

(1)  For  stationary  or  windowed  signals  find  the 
correlation  vector  r^n\  For  nonstationary 
systems  or  infinite  unwindowed  data  find  the 
correlation  matrix  R^n\ 

(2)  Find  v/n^.  For  the  first  stage  of  analysis  the 
initial  conditions  of  (4-16)  is  used.  For 
later  stages  (4-13)  gives  the  recursive  solution 
for  w(n). 

(3)  The  numerator  of  (4-30)  is  found  using  the  weight 
vector  and  either  a  correlation  vector  r^n^  or 
the  correlation  matrix 
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Figure  4.2  Flow  Chart  for  mth  Order  Analysis 


L  35 


. 


(4)  Find  P(n)  from  (4-28).  For  n=0,  PCo)  =  Ryy(0); 
otherwise  use  the  recursion  of  (4-28). 

(5)  ^tn+l)  is  now  found  from  (4-30). 

(6)  If  higher  order  analysis  is  desired,  increment 
n  and  return  to  (2). 

Figure  4.3  gives  a  comparison  of  the  number  of 
multiplications  required  to  calculate  reflection  factors  for 
a  seven  stage  lattice  that  uses  256  data  points,  and  assumes 
a  stationary  system  with  the  data  windowed  to  ensure  a 
Toeplitz  structure. 


CONDITIONS:  7  stages 

256  data  points 

Data  windowed  for  Toeplitz  and  symmetric 
correlation  matrix 

TIME  DOMAIN  LATTICE: 

#  multiplications  =  1+  (£+2)n  +  Mn-1)2  = 

p(o)  error  signal 

update 

256  +  (256+2)7  +  256(6)2  =  5134 

For  Correlation  Lattice  ,  „  0  _ 

m-1  m  /  ., v i  m 

#  multiplications  =  l  +  Z  (n+1)  +  Z  — = -  +  m£.+  Z  (n+1) 

n=0  n=l  n=l 

R  (0)  w(n)  =  w(n)  R  (n)  k(n) 

yy  -  yy 

=  256  +  28  +  140  +  1792  +  34  =  2250 


Figure  4.3 
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Multiplications  Required  for 
Implementation 
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B.  THE  TWO  CHANNEL  LATTICE 

The  two  channel  lattice  developed  in  II, C  can  similarly 
be  redefined  in  terms  of  correlation  weights,  the  details  of 
which  are  presented  in  Appendix  A.  The  following  is  a 
summary  of  the  steps  and  equations  necessary  to  implement 
the  algorithm. 

To  implement  the  two  channel  algorithm  it  is  necessary 
to  calculate  eight  weighting  vectors,  two  for  each  of  the 
four  forward  and  backward  error  signals.  These  can  be  found 
using  the  initial  conditions 
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and  the  recursions  for  higher  orders 
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and 
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The  vectors  of  (4-33)  and  (4-34)  each  have  two  elements 
which  are  also  vectors,  and  the  matrix  multiplication  is 
executed  as  a  (2x2)  matrix  times  a  (2x1)  vector. 

The  (2x2)  matrix  A* 
vectors  and  a,  correlation  matrix  Rv 
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can  be  found  using  these  weight 
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The  four  elements  of  the  matrix  are  defined  by  the 

relations 
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The  P  and  P  matrices  are  found  using  the  recursions 
previously  cited,  namely 


P(0)  =  p(0)  s 
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The  lattice  parameters  are  found  using  the  A^n^ ,  p^n^ 
and  p^n^  matrices  from 


(n+1)  p(n)  1  A(n)T 
K  =  P  A 


(4-44) 
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A  flow  chart  of  this  implementation  is  presented  as 
Figure  4.4. 


[i] 


Figure  4.4  The  m-order  Two  Channel  Implementation 
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Correlation  Lattice  Simulations 

To  verify  the  algorithms  presented  in  the  forgoing  a 
FORTRAN  program  was  written.  This  program  implements  the 
flowchart  of  Figure  4.4.  The  correlation  matrices  generated 
by  this  program  are  both  Toeplitz  and  symmetric.  The 
examples  here  are  the  same  two  systems  used  in  Chapter  III.E. 
and  represent  a  run  which  is  typical  of  the  ensemble  of  tests 
that  were  made .  The  distribution  over  the  ensemble  was  minor 
since  measurement  noise  was  not  considered.  The  fourth  order 
model  of  System  1  is  shown  with  Gaussian  elimination  using 
256  points,  in  Figure  4.5.  Figure  4.5  shows  the  results 
using  4096  points.  Figure  4.5  shows  the  same  comparison  of 
System  2  using  128  data  points,  and  Figure  4.8  is  System  2 
modeled  with  1024  points.  In  almost  all  cases  the  results  of 
the  correlation  method  are  essentially  identical  to  the 
results  of  the  Gaussian  elimination  method.  This  should  not 
be  unexpected.  Both  matrix  manipulation  and  correlation 
lattice  methods  use  identical  data,  the  auto  and  cross 
correlations  of  the  input  data  signals.  Both  methods  window 
the  data  identically.  Both  methods  are  solving  the  same 
normal  equations,  though  through  differing  algorithms. 

The  number  of  multiplications  saved  by  the  correla¬ 
tion  lattice,  when  compared  to  the  conventional  lattice 
implementation,  is  approximately  8p,  where  p  is  the  number 
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of  data  points.  This  is  because  the  error  signals  are  not 
multiplied  by  the  lattice  parameter  matrices  and 

in  the  correlation  implementation. 
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V.  SUMMARY,  CONCLUSIONS  AND  OPEN  QUESTIONS 

The  objective  of  this  dissertation  has  been  to  develop 
efficient  methods  of  implementing  the  AR  lattice.  Two 
particular  applications  for  system  identification  and 
modeling  have  provided  the  motivation  for  this  research. 

The  first  is  system  fault  identification.  Computationally 
efficient  algorithms  must  be  implemented  if  the  fault 
analysis  is  expected  to  be  done  on  a  real  time  system. 
Adaptive  lattice  structures  require  only  moderate  computation 
for  each  data  point  and  have  often  been  implemented  in  real 
time.  The  block  structures,  such  as  the  conventional  time 
domain  or  frequency  domain  lattices,  require  greater 
numbers  of  calculations  in  bursts  for  each  block  of  data. 
These  structures  could  be  implemented  using  both  parallel 
and  pipeline  procedures:  parallel  calculation  of  the 
correlations  of  each  channel  and  pipelined  calculation  of 
the  stages  of  the  model. 

The  correlation  lattice  structure  also  encourages 
parallel  processing  implementations.  In  this  structure 
calculation  of  the  correlation  matrix  requires  the  greatest 
computational  effort.  These  correlation  estimates  can  be 
performed  in  parallel  as  the  data  becomes  available. 

The  second  application  considered  is  the  modeling  of 
nonlinear  systems.  For  this  application  the  foregoing 
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issues  are  of  even  greater  importance.  Hybrid  signals, 
products  and  cross  products  of  the  nonlinear  system  input 
and  output,  are  generated  and  the  resultant  signals  applied 
as  the  input  of  a  multichannel  lattice.  In  this  case  the 
computational  and  storage  savings  of  the  correlation  lattice 
become  of  major  significance  and  allow  the  implementation  of 
more  complex  models.  It  has  been  found  that  the  two  channel 
lattice  when  applied  to  ARMA  modeling  can  be  simplified 
because  of  the  relationships  that  exist  between  the  input 
and  the  output  of  the  system  under  test.  The  lattice  can 
be  simplified  further  when  the  system  driving  signal  is  white. 

Some  of  the  constraints  found  are  given  physical  inter¬ 
pretations  that  relate  to  the  system  under  analysis.  There 
has  been  no  such  interpretation  for  the  relationship  of 
the  lattice  parameters  by  a  factor  of  the  DC  gain  of  the 
system  being  modeled.  It  is  felt  that  this  simple  relationship 
should  have  some  physical  interpretation  related  to  the 
system  under  test. 

The  lattice  was  redefined  using  the  frequency  domain 
representation  of  the  input  data.  This  was  found  to  be 
useful  because  operations  in  the  frequency  domain  allowed 
for  normalization  of  the  input  signals  to  whiten  their 
spectra.  This  whitened  spectra  provides  a  standard  input 
that  allows  lattice  parameters  to  be  used  for  identification 
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and  ensures  the  previously  mentioned  relationships  of  the 
lattice  parameters . 

There  are  two  normalization  procedures  proposed.  The 
first  requires  the  averaging  of  ensembles  of  the  inputs  to 
obtain  an  expectation  of  the  energy  spectrum  of  the  input. 
The  second  method  uses  the  spectrum  of  a  single  sample 
spectrum  as  a  normalization  function.  Either  of  these 
methods  will  fail  if  the  energy  at  any  single  spectral 
harmonic  is  zero . 

It  is  not  understood  how  either  of  these  methods  affects 
the  accuracy  of  the  resultant  model.  Experiments  comparing 
normalization  methods,  with  and  without  spectral  smoothing 
introduced  by  windowing,  have  been  inconclusive.  While  some 
examples  have  shown  normalization  can  improve  model  accuracy 
over  an  unnormalized  analysis,  other  examples  show  the 
opposite.  It  is  postulated  that  the  dynamic  range  of  the 
input  spectrum  and  normalization  spectrum  are  of  major 
significance,  but  there  has  been  no  analysis  of  the 
procedures  involved. 

It  is  also  found  that  this  frequency  domain  lattice  can 
be  used  to  synthesize  a  digital  filter.  The  desired 
frequency  domain  spectrum  is  provided  as  the  input  of  the 
lattice,  and  recursive  in  order  lattice  algorithm  provides 
an  LMS  lattice  model  which  meets  the  desired  specification. 


Lastly  it  is  found  that  the  lattice  can  be  reformulated 
in  terms  of  correlation  weights  of  the  input  signals.  This 
is  efficient  because  it  becomes  unnecessary  to  generate 
error  signals  within  the  lattice.  This  results  in  both 
computational  and  storage  efficiencies  when  compared  to 
conventional  lattice  implementations.  Though  developed  for 
the  two  channel  lattice,  the  procedure  can  be  extended  to 
higher  dimensioned  lattices  where  the  efficiencies  will  be 
of  even  greater  significance. 

Inherent  in  the  derivation  of  Levinson's  algorithm  is 
the  requirement  that  the  input  signals  correlation  matrix  be 
symmetric  Toeplitz,  requiring  that  the  signals  be  stationary 
and  correlations  estimates  are  from  windowed  data.  When  the 
lattice  structure  is  derived  from  Levinson's  algorithm  the 
requirement  of  windowing  the  data  is  displaced  by  the 
lattice  structure,  which  forces  upon  the  data  signals  a 
different  but  nevertheless  present  window  function.  The 
symmetric  Toeplitz  nature  of  the  correlation  matrix  is  not 
an  issue  with  the  lattice  structure. 

When  the  lattice  is  reformulated  in  terms  of  correlation 
weights,  it  is  again  necessary  to  generate  a  correlation 
matrix.  It  does  not  appear  in  the  derivation  that  there  is 
a  requirement  for  a  symmetric  Toeplitz  matrix  but  it  is  also 
not  understood  how  a  nonsymmetric  Toeplitz  correlation 
matrix  will  affect  the  accuracy  of  the  solution.  This 
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formulation  can  be  expected  to  generate  errors  in  the 
solution,  and  it  is  not  clear  how  they  will  propagate 
through  the  correlation  lattice  algorithm.  This  is  an  issue 
that  deserve  future  attention. 


APPENDIX  A 


The  two  channel  lattice  structure  of  Figure  2.8  can  be 
reconfigured  for  a  more  efficient  implementation  using 
correlation  weights  rather  than  error  signals.  The  lattice 
equations  earlier  defined  in  II. C  are  repeated  here  for 
convenience . 

( A-l ) 

(A-2 ) 

e(n+1)(k)  =  e(n)(k)  -  K(n+1)  e(n)(k-l)  (A-3) 

e(n+1)(k)  =  e<n)(k-l)  -  K(n+1)  e(n)(k)  (A-4) 

and  the  reflection  factor  solutions 

P<0)  =  P<0)  =  s{e<0)(k)e(0,T(k»  = 
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p(n+l)  __  p C n ) ^ ^ ( n ) ^ ( n ) ^ 


(A-7 ) 


A(n)  =  e{e(n)(k)e(n)(k-l)}  (A-8) 
K(n+D  =  p(n)  (A-9 ) 
K(n+1)  =  p(n)  1A(n)  (A-10 


A  channel  error  signal  e^n^(k)  or  e^n^(k)  is  a  linear 

X  X 

combination  of  the  present  value  and  the  past  n  values  of  y 
and  u.  Using  this,  e^n^(k)  may  be  written 

>><*[  = 

e(n)(k) 
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Vectors  y^n\k)  u^n\k)  and  w  are  each  (n+1)  terms  long. 

—  — xz 

The  result  of  the  multiplication  of  (A-12)  is  a  2  element 


column  vector. 


Using  this  same  notation  the  backward  delayed 
vector  e  (k-1)  may  be  written 
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(A-17) 


where 


(A-18) 


In  the  single  channel  case  the  backward  weighting  vector 
was  equal  to  the  forward  weighting  vector  inverted  in  place. 
This  arose  because  forward  and  backward  reflection  factors 
are  equal  in  the  single  channel  lattice.  In  the  two  channel 
lattice  is  not  in  general  equal  to  so  there  is  no 

similar  relationship  between  the  forward  and  backward 
weighting  vectors. 
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This  provides  a  recursive  solution  for  the  forward  weight 
vectors . 

In  illustration,  e^(k)  is  defined  by  (A-l)  as 


(A-21) 


=  CyCJOyCk- 
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This  result  can  be  substituted  into  (A«ll) 
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(A- 2 4 ) 


Inspection  of  Figure  2.8  shows  that  (A-24)  is  indeed  the 
equation  of  the  first  stage  error  signal. 

An  analogous  procedure  can  be  allowed  to  develop  the 
recursive  solution  for  the  backward  weights.  Substituting 
(A-12)  and  (A-17)  into  (A-4)  yields 


e(n+1)(k)  =  [£(n+1)T(k)j uCn+1)T(k)] 
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With  the  weight  vectors  for  all  stages  established,  the 
reflection  factors  may  be  found  in  terms  of  these  weights. 
Substituting  (A-12)  and  (A-17)  into  (A-8)  yields 


The  center  matrix  is  dimensioned  (2x2).  Its  four  elements 
are  each  a  submatrix  of  dimension  C2Cn+2)x2(n+2) ) ,  Each 
submatrix  is  premultiplied  by  a  2(n+2)  element  row  vector 
and  postmultiplied  by  a  2Cn+2)  element  column  vector.  The 
result  of  these  multiplications  will  be  a  scalar  value  for 
each  of  the  four  elements  of 

The  order  of  these  vector  multiplications  may  be 
inverted  yielding 
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Extracting  An  and  rewritting 
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is  a  (n+2)x(n+2)  (cross)  correlation  array. 
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By  defining  a  2 Cn+2 )x2 (n+2 )  matrix  R(ni  as 
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Implementation  is  discussed  in  Chapter  IV. 3. 
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